Host–pathogen dynamics in longitudinal clinical specimens from patients with COVID-19

Rapid dissemination of SARS-CoV-2 sequencing data to public repositories has enabled widespread study of viral genomes, but studies of longitudinal specimens from infected persons are relatively limited. Analysis of longitudinal specimens enables understanding of how host immune pressures drive viral evolution in vivo. Here we performed sequencing of 49 longitudinal SARS-CoV-2-positive samples from 20 patients in Washington State collected between March and September of 2020. Viral loads declined over time with an average increase in RT-QPCR cycle threshold of 0.87 per day. We found that there was negligible change in SARS-CoV-2 consensus sequences over time, but identified a number of nonsynonymous variants at low frequencies across the genome. We observed enrichment for a relatively small number of these variants, all of which are now seen in consensus genomes across the globe at low prevalence. In one patient, we saw rapid emergence of various low-level deletion variants at the N-terminal domain of the spike glycoprotein, some of which have previously been shown to be associated with reduced neutralization potency from sera. In a subset of samples that were sequenced using metagenomic methods, differential gene expression analysis showed a downregulation of cytoskeletal genes that was consistent with a loss of ciliated epithelium during infection and recovery. We also identified co-occurrence of bacterial species in samples from multiple hospitalized individuals. These results demonstrate that the intrahost genetic composition of SARS-CoV-2 is dynamic during the course of COVID-19, and highlight the need for continued surveillance and deep sequencing of minor variants.

www.nature.com/scientificreports/ co-occurrences, which may be associated with recovery 12,13 . We found negligible change in viral consensus sequences over time, but detectable changes in variant allele frequencies that are only weakly predictive of future consensus changes across the globe. We further observed rapid emergence of deletion variants in the N-terminus domain of the spike glycoprotein in one patient, potentially suggesting within-host SARS-CoV-2 evasion of NTD-directed antibodies. Taken together our results support the limited emergence and fixation of escape variants during a typical infection, and also highlight the need to monitor minor variants due to their potential impact on vaccine and therapeutic efficacy.

Viral load dynamics in longitudinal samples from SARS-CoV-2 infected individuals. Resid-
ual clinical specimens were obtained from the University of Washington (UW) Virology Lab after testing for SARS-CoV-2 14 . By reviewing our laboratory information system, we identified 20 individuals who had two or more positive or inconclusive samples collected between March and September 2020 (Table 1, Supplementary  Table S1). Inconclusive samples had one of two PCR targets detected. This commonly occurs in samples with small amounts of viral DNA, and thus these were considered presumptive positive for SARS-CoV-2, albeit with low viral load. A majority of samples came from inpatients who received care within the UW Medicine system, which includes UW Medical Center, Harborview Medical Center, and Northwest Hospital. Samples came from individuals with a mean age of 70 (range 42-99), many with severe disease given the availability of multiple samples from these patients. Consistent with other reports 15 , we observed that viral load declined over time in most patients with two or more positive or inconclusive samples with an average increase in RT-QPCR cycle threshold (Ct) of 0.87 per day (Fig. 1, Supplementary Fig. S1). A total of 49 samples (47 nasopharyngeal and 2 oropharyngeal swabs) were sequenced with sufficient reads to be included in this study ( Supplementary Fig. S2A). Ct values for these samples ranged between 14.6 and 36.6. The length of time between collection dates for sequenced samples from the same individual ranged from 0 to 22 days. Samples collected on the same date were sequenced for three individuals. We sequenced nasopharyngeal and oropharyngeal samples from P004 collected at the time of autopsy, two samples from P006 collected at the same time during a hospital admission, and samples collected 9 h apart from P012 during an emergency room visit.
Negligible change in consensus sequences over time. We obtained full-length viral genome sequences with less than 2% unknown bases (Ns) for two or more time points in 14 out of 20 individuals, Low frequency variants detected across the genome. We analyzed intrahost viral genetic variation by examining all sites with > 100 × locus depth, masking known problematic sites and filtering for high-confidence variants (see Methods). We examined sites in 47 samples from 20 different patients and found a total of 103 unique non-synonymous variants relative to the Wuhan-Hu-1 (NC_045512.2) reference genome present at frequencies between 5 and 95% ( Fig. 2A). nsp12 had the highest number of variant sites (26, 25.24%), followed by nsp3 (16,15.533%) and the spike glycoprotein (13, 12.62%). Even when adjusted for gene length, mutations were most prevalent in nsp12, at 0.009 variant sites per nucleotide. For samples that were sequenced multiple times (n = 10), variant frequencies were reproducible across replicates, particularly among samples sequenced with the same library preparation method (Fig. 2B). Intra-host mutations present in multiple longitudinal samples of the same patient were mostly fixed and highly clonal, with some low-level variants ( Supplementary Fig. S2B).
Of the seven most commonly observed variants in our dataset (Table 3), the three most frequent define the Washington state outbreak clade 16,17 and the rest of the variants are clade-defining mutations in Nextstrain clades 20A and 20C. Nine out of 20 patients had the spike protein mutation D614G (A23403G), which has been associated with increased transmissibility and higher viral loads 18,19 . While this variant was rare at the beginning of the pandemic, it reached near fixation in the global SARS-CoV-2 population by June 2020 20 . This rapid rise in prevalence is reflected in our data, as this D614G mutation is present in all three patients with samples collected during or after June 2020. In 10 out of the 11 patients with the 614D variant, no alternate alleles were detected at this position. In the second sample from P007, 614G was detected with a variant allele frequency of 6.1%, but the read depth at this locus (82X) was insufficient to reach our QC standards.
Variants exhibiting intra-host evolution are limited in prediction of future global consensus changes and highlight SARS-CoV-2 antibody evasion. We further examined variants that underwent a maximum allele frequency change of ≥ 20% across timepoints within each patient. The derived alleles for all 25 non-synonymous amino acid changes meeting this criteria in our dataset were also observed among consensus genomes deposited in GISAID 21 by April 2021 (range: 1-2171, mean: 467.8, median: 86.5). Only 8 variants exhibited a maximum allele frequency change of ≥ 40%. The derived alleles for these variants were present in very few GISAID consensus sequences (range: 1-1751), which at the time of analysis represented a mere 0.0001-0.17% of all sequences deposited in GISAID (Fig. 3A). Though the number of consensus sequences with these derived alleles was relatively low, these sequences were diverse with respect to collection date and geographic origin (Fig. 3B).
Three variants with a ≥ 20% within-host change in allele frequency localized to the spike glycoprotein (Fig. 4A). None of these had an allele frequency change of ≥ 50%. In patient P016, a S:143-145 6-nucleotide deletion was observed at an allele frequency of 20.7% in patient sample 2, but was not observed in samples collected www.nature.com/scientificreports/ 3 days prior and 3 days later. This patient was immunocompetent and not receiving any COVID-19 treatment. Interestingly, numerous other deletions arose at low frequencies in this patient, with the largest number present at day 0 (Fig. 4B). The most prevalent deletion variant in the day 0 sample (collected 16 days after initial symptom onset) was S:∆141-144 at 1.87% allele frequency. This deletion was the second most common variant in the day 3 sample at 1.2% allele frequency, but was not observed at all in the day 6 sample. Deletions in this region, including S:∆141-144, have previously been observed in chronically infected immunocompromised patients, and some are associated with escape from NTD-specific neutralizing antibodies or polyclonal sera [8][9][10][11]22 . Notably, a S:∆142-144 deletion is a hallmark of the currently circulating Omicron variant.

Longitudinal RNAseq analysis illustrates loss of ciliated epithelium during infection. For sam-
ples that were sequenced metagenomically, we pseudo-aligned reads to the human transcriptome to perform differential expression analysis comparing initial (t = 0) timepoints to later timepoints. Samples with more than 900,000 pseudo-aligned reads (n = 7 initial, 3 later timepoints) were included in the analysis to determine vari- Table 2. Consensus sequence analysis of SARS-CoV-2 in longitudinal specimens. All patients with less than 2% unknown bases (Ns) are included. The last column indicates nucleotide differences compared to the first sample collected for each respective patient. One asterisk (*) indicates symptoms were present at first time point but exact date of symptom onset is unknown. Two asterisks (**) indicate days since first sample.  www.nature.com/scientificreports/ ation in host gene expression over time. We observed a dramatic downregulation of several cytoskeletal genes, particularly dynein heavy chain (DNAH 2, 3,5,6,7,9,10,11,12), as well as WDRs, MAP1A, and others (Fig. 5A, Supplementary Table S3). Gene Ontology analysis (Fig. 5B) confirmed that downregulated genes are involved in biological processes associated with microtubule-based motility. This is consistent with the death of ciliated epithelial cells, which are enriched for transcripts encoding microtubule transport machinery 23 , following SARS-CoV-2 infection. We observed upregulation of some actin cytoskeleton-related transcripts like VAV1, VASP, and RhoF 24 . We also observed downregulation of several interferon-stimulated genes, but these did not reach statistical significance in this small sample set.
Metagenomic analysis shows high levels of clinically relevant bacteria in three samples. We used a previously described metagenomic pipeline (CLOMP 25 ) to perform taxonomic assignment of preprocessed reads. We excluded samples that had fewer than 10,000 reads after trimming and any samples that underwent enrichment via probe-capture or amplicon-sequencing for SARS-CoV-2. A total of 24 samples from 11 patients were included in further analysis ( Supplementary Fig. S3, Supplementary Fig. S4). No viruses aside from SARS-related coronaviruses met the required cutoffs to be classified as co-infections (see "Methods"). Multiple samples had detectable numbers of bacterial reads, in particular Staphylococcus aureus and Moraxella catarrhalis. Samples from two patients (P004, P005) had high detectable levels (> 100 RPM) of both bacteria at one or more time point(s) (Supplementary Fig. S3B). In patient P006, who had paired nasopharyngeal (NP) and oropharyngeal (OP) swabs collected at the same time point, we detected Capnocytophaga gingivalis, Capnocytophaga leadbetteri, and Streptococcus parasanguinis in the OP swab at 21,960, 9475, and 4476 RPM, respectively. All three species of bacteria commonly colonize the oropharynx. In contrast, in the NP swab, the predominant species of bacteria was a common skin colonizer, Cutibacterium acnes, at 824 RPM. In P010, we found a large number of reads corresponding to Corynebacterium spp. at one time point. Upon review of medical records, we found no mention of bacterial co-infections or of positive bacterial cultures in the charts for P004 or P006. P005 had a nares culture that grew methicillin-resistant Staphylococcus aureus three months prior to SARS-CoV-2 infection.

Discussion
In this study, we performed high-throughput sequencing of longitudinal clinical specimens that were positive for SARS-CoV-2 by RT-QPCR. Most samples were sequenced using a metagenomic approach, which enabled us to simultaneously derive information about viral evolution, host transcription, and the presence of other organisms within patient samples. We showed that although the viral consensus sequence remains largely unchanged over the course of infection, there is a relative abundance of genome-wide low-frequency variants. Similar to other studies, we saw a wide range in the number of variants detected across samples 26 and distribution of variants across the genome, though some positions appeared to be more prone to variation 4,26 . Studies of SARS-CoV-2 and other respiratory viruses 2,3,8,27 have demonstrated the transmission of minor variants and the role of these population bottlenecks on viral evolution, underscoring the importance of studying within-host viral variation. All variants demonstrating significant longitudinal evolution in our sample set collected March-September 2020 have been observed in consensus sequences from around the globe 21 , albeit at relatively low prevalence.
In one patient, we observed rapid turnover of multiple deletion variants in the N-terminal domain of the spike glycoprotein, which has previously been seen in persistent infection in immunocompromised individuals and has been associated with viral escape of neutralizing antibodies 8,10,22,28 . Deletions in the NTD are of particular significance due to their presence in currently circulating lineages of concern. Here we show the emergence of a deletion in this genomic region in an immunocompetent background. It is unclear if the absence of this mutation at day 6 is due to successful clearance of the NTD variant or lack of detection of the minor allele associated with lower copy numbers. In addition, while some of these low frequency deletion alleles have been previously shown to arise independently in different patients in response to similar selection pressures 8 , the presence of multiple low frequency deletion alleles within the same patient may be the product of parallel within-host microevolutionary processes suggestive of some selective advantage, particularly with its similarity to the S:∆142-144 deletion in the currently circulating Omicron variant. Notably, we did not find any evolving www.nature.com/scientificreports/ variants selected for in the RBD, the main target for neutralizing activity of human plasma 29,30 . Additionally, we found two highly dynamic mutations within the Ubl2 and PL2Pro domains of the multifunctional nsp3 protein 31 and three mutations within the nsp12 polymerase. While we and others have previously linked a specific nsp12 mutation to remdesivir resistance 32 , each of these mutations has no phenotypic associations to date and awaits further biochemical and virological characterization. Individual host factors, such as the immune response and respiratory tract microbiome, may play an important role in viral persistence. In particular, because SARS-CoV-2 infection is slow to resolve, the adaptive immune response could drive within-host viral evolution as variants that can escape T-cell and antibody responses develop. Although we were underpowered to see specific evidence of an adaptive immune response being mounted against SARS-CoV-2 in the nasopharynx, detailed studies evaluating antibody and T-cell receptor repertoire changes throughout the course of infection could shed light on the role of immune pressure in the development of minor variants. Similarly, the relationship between bacterial colonization of the nasopharynx  . Variants that exhibit intra-host evolution in the spike protein across all patients. (A) All nonsynonymous variants located in the spike protein with a ≥ 20% change in allele frequency among timepoints for any patient. (B) Enumeration of deletions that arose between residues 138-149 of the spike protein in P016 at ≥ 1% relative frequency reveals a rapidly changing complement of low frequency alleles present over a 6-day period. The reference nucleotide sequence (NC_045512) is located at the top of the sequence alignment. Above the reference is the corresponding amino acid sequence with associated residue numbers. Alleles that match the reference are in gray, and deletions are shown in black. To the right of the sequence alignment is a bar graph showing the square root of the relative frequency of each variant, for visualization purposes, labeled with the allele frequency in percentage and read depth in parentheses. www.nature.com/scientificreports/ and the development or suppression of inflammation in response to SARS-CoV-2 infection remains poorly understood.
As viral load decreases during recovery, it becomes more challenging to recover viral genomes. As a result, one of the limitations of our study is the variability in sequencing depth across samples and the difficulty in ensuring similar sequencing depth for samples from different time points. We used an amplicon sequencingbased approach described previously 33 to obtain near full-length genomes from low viral load samples (up to Ct values of 36). We also used multiple library preparations and performed re-sequencing to ensure the accuracy of variant calls.
Taken together, our results suggest that low frequency genomic variants emerge in immunocompetent individuals, but that these variants are unlikely to reach fixation. Given the emergence of rapidly spreading variants of concern over the past several months, the limited intra-host evolution observed in our dataset highlights the critical impact that a select few individual intra-host evolutionary events may have on the course of the global pandemic and the need for continual genomic surveillance.

Methods
Sample collection and clinical testing for SARS-CoV-2. Specimens were obtained as part of clinical testing for SARS-CoV-2 ordered by local healthcare providers or collected at drive-through testing sites. RNA was extracted and the presence of SARS-CoV-2 was detected by RT-QPCR as previously described using either the emergency use-authorized UW CDC-based laboratory-developed test, Hologic Panther Fusion or Roche cobas SARS-CoV-2 tests 34,35 . Supplemental Table 1 contains clinically relevant details of these specimens.
Chart review and ethics approval. We received approval from the University of Washington Institutional Review Board (UW IRB) to use residual clinical specimens for sequencing and to review clinical records of patients who received care within the UW network. We also received a waiver of informed consent from the UW IRB for the use of residual clinical specimens and retrospective chart review to perform this work. Information obtained from medical records included sex, age, comorbidities, medication, hospital or critical care admission, and discharge status. All methods were carried out in accordance with relevant guidelines and regulations. 18 of the 20 total patients had one or more known comorbidities (Supplementary Table 1).
Sequencing and bioinformatic analysis. Sequencing was attempted on all samples with a positive RT-QPCR assay result that had a Ct ≤ 36 using one of three methods available at the time: (1) a shotgun metagenomic approach using Illumina Nextera XT described previously 36 for samples with a Ct less than 24; (2) an oligo-nucleotide probe capture-based approach from IDT (xGen NGS hybridization capture, Integrated DNA Technologies) similar to previous work 37 for samples with Ct between 24 and 28, or (3) using Swift Biosciences' Normalase Amplicon Panel library preparation workflow 38 for samples with Ct between 28 and 36. Libraries were sequenced on Illumina MiSeq, NextSeq, or NovaSeq instruments using 300, 150, 100, or 75 bp reads. Consensus sequences were assembled using TAYLOR 38 , a custom bioinformatics pipeline (https:// github. com/ greni nger-lab/ covid_ swift_ pipel ine) with or without an additional primer clipping step depending on library preparation method.
Consensus sequences from each individual were aligned with the reference sequence NC_045512 using MAFFT v7 39 . Clade assignments were generated using Pangolin (http:// github. com/ cov-linea ges/ pango lin) and Nextstrain 40 in December 2020. Consensus sequences with < 5% Ns across the length of the genome were considered for further analysis.
Variants were also called with TAYLOR from aligned reads. Variants leading to coding changes with a sequencing depth of > 100 and an allele frequency > 0.01 were subjected to further analysis. For variants that had an allele frequency of < 0.4, we included an additional quality-control step to only include mutations that were present across multiple samples. We also excluded mutations in the first 100 and last 50 bases, as well as variants determined to be due to sequencing error. Most samples were re-prepped and sequenced multiple times to ensure accuracy of variant calls. Variants at positions 6700, 11081-83, 19989, and 29056 were observed in a large number of samples but were determined to be the result of homopolymer sequencing error and were excluded.

RNAseq analysis.
For samples that were prepared with no target enrichment step (metagenomic sequencing), reads were adapter and quality trimmed with Trimmomatic v0.39 41 using the call "leading 3 trailing 3 slidingwindow:4:15 minlen 20", then pseudoaligned to the hg38-derived human transcriptome using Kallisto v0.46 42 . Only samples with more than 900,000 reads pseudo-aligned to the human genome were used for analysis. Differential expression analysis using the Wald test was performed using DEseq2 43 and deemed significant at a Benjamini-Hochberg adjusted p value < 0.1. Statistical enrichment of Gene Ontology Biological Processes was performed on all significant genes using the R package clusterProfiler 44 . Raw counts were submitted to the Gene Expression Omnibus, accession number GSE173310.
Metagenomic analysis. Raw FASTQ files were analyzed using CLOMP v0.1.4 (https:// github. com/ FredH utch/ CLOMP) as previously described 25 . Samples with more than 10 million reads were randomly down-sampled to 10 million reads before analysis using the "sample" command in seqtk (https:// github. com/ lh3/ seqtk). The pipeline output was visualized using the Pavian metagenomic explorer 45 , and reads per million (RPM) calculations were done using a custom R script. Results were filtered to highlight RPM counts for a shortlist of clinically relevant taxa (Supplementary Table S4). Samples were determined to be positive if the species level RPM was at least 30 for viruses, and 100 for bacteria.


Comparative mitochondrial genome analyses reveal conserved gene arrangement but massive expansion/contraction in two closely related Exserohilum pathogens

Graphical abstract

Introduction
Exserohilum is an important plant-pathogenic genus belonging to the order Pleosporales (Dothideomycetes, Ascomycota), and until now, a total of 40 Exserohilum species have been recorded (https://www.indexfungorum.org). This taxonomic group was originally described as Helminthosporium, which was then subdivided into Bipolaris, Curvularia and Exserohilum based on the traditional taxonomical criteria especially by their hilum morphology [1][2][3]. The classification result was also supported by multigene phylogenetic analyses based on nine nuclear loci including ITS, LSU, act, tub2, cam, gapdh, his, tef1 and rpb2 [3]. Although known as Setosphaeria (teleomorphic stage) as well, the scientific name of anamorph Exserohilum was generally recommended for use according to the International Code of Nomenclature for algae, fungi and plants [3,4].
As the type species of Exserohilum, E. turcicum, a filamentous hemibiotrophic fungus, is notorious for causing northern corn leaf blight (NCLB) [3,5]. Of late years, this disease has become increasingly problematic worldwide with the extensive planting of susceptible varieties of maize in many countries [6,7]. The most common symptom of NCLB is necrotic green-gray lesions on leaves with a length range from 2 to 14 cm, which could lead to 10-40% grain yield losses in severe cases [8]. Besides is also an important blight agent of some Sorghum spp. [9]. By contrast, another noteworthy species of this genus, E. rostratum, has a more complex life style. It was first described as the pathogen of thrush grass in 1923 in the USA [10], and then reported to infect more than 30 plant species belonging to 28 genera of 11 families [11], including gramineous crops like corn and rice [12], woody plants like rubber tree [13], and some tropical cash crops like sugarcane [14]. Recently, E. rostratum seemed to be better known as an opportunistic pathogen found in humans, causing keratitis, meningitis, subdermal fungal disease, etc. [15][16][17]. Between the above two closely related Exserohilum species with fairly different target hosts, comparative analysis was rarely performed. Even if their nuclear genomes have been already available, it is not clear how E. turcicum and E. rostratum differ in their virulence strategies and molecular evolution mechanisms [15,18].
Mitochondria are a kind of double membrane organelles wellknown for their major role in energy supply, and therefore commonly referred to as the power factories of eukaryotes. In addition, they can also participate in the ion homeostasis, intermediate metabolism, cell senescence and apoptosis [19]. In recent years, an increasing body of evidence suggested that mitochondria were derived from the symbiotic combination of a-proteobacteria and eukaryotic cells [20,21]. And their evolutionary processes varied among the main eukaryote groups including animals, plants and fungi [22], of which, however, fungal mitochondria shared the least spotlight. The mitochondrial genomes (mitogenomes) of different fungal species showed significant differences in genome size, gene arrangement, ratio of repetitive and intergenic regions, and number of introns and open reading frames (ORFs) [23,24]. Among all published fungal mitogenomes, Hanseniaspora uvarum has the smallest one that is only 18.84 kb in size [25], while the mitogenome of Morchella crassipes is the largest, approximately 531.19 kb in size [26]. In spite of the varying sequence lengths, fungal mitogenomes usually contain a set of conserved protein-coding genes (PCGs) including seven NADH deaminase complex subunits (nad1-6 and nad4L), three cytochrome c oxidase subunits (cox1-3), one cytochrome b gene (cob) and three ATP synthase subunits (atp6, atp8 and atp9), two ribosomal RNA (rRNA) genes, and 22-26 transfer RNA (tRNA) genes [27,28]. Studies have found that expansion and contraction of the mitogenomes in different fungi were mainly due to the number and length of introns, repeat rate and introduction of new genes through horizontal transfer [29,30]. There are two types of introns in fungal mitogenomes, namely group I and group II, in which homing endonuclease genes (HEGs) were often detected and could produce LAGLIDADG/GIY-YIG motif endonucleases [31,32]. With the development of the nextgeneration sequencing (NGS) technology and continuous reduction of sequencing costs, more than 12,500 complete mitogenome sequences of eukaryotes have been released to date in the NCBI database (https://www.ncbi.nlm.nih.gov/genome/browse/#!/organelles/). Nonetheless, for most fungal taxa including the genus Exserohilum, there has been no research on their mitochondrial genomes, which will limit our knowledge of the ''second genome" in some significant fungal lineages.
In the current study, mitogenomes of the two Exserohilum pathogens, E. turcicum and E. rostratum, were sequenced, assembled, annotated, and compared with some other Pleosporales species. Our main objectives are as follows: 1) to depict the mitogenome contents, structures and organizations of the two Exserohilum species; 2) to perform comparative mitogenomic analysis between the different Exserohilum species and reveal their variations and conservations; 3) to reveal the intron dynamic changes of cox1 genes in 13 Pleosporales species including the two Exserohilum fungi; 4) to clarify the phylogenetic status of Exserohilum in the phylum Ascomycota based on combined mitochondrial gene sets. This study served as the first investigation on Exserohilum mitogenomes, which will lay the foundation for further understanding the genetic evolution, species differentiation, and ecological adaptation of Exserohilum species.

Fungal isolates, DNA extraction and genome sequencing
With the method of tissue isolation [11], both E. turcicum and E. rostratum strains (ZM10601 and ZM170581) were obtained from corn leaves in the city of Zhengzhou (34°48 0 N, 113°39 0 E), Henan province, China. Their hyphal tips were transferred onto PDA plates for purification, and after sufficient growth, these culture plates were kept at 4 ℃ for short-term storage. For long-term preservation, however, the fungal colonies needed to be rinsed off twice with sterile water for collecting conidia, which were then stored in 15% glycerol at À 80 ℃ in the fungal collection of Henan Agricultural University. Species identification was based on their morphological characters and sequence analysis of the ITS, LSU, gapdh, tef1 and rpb2 genes. These sequences have been deposited in GenBank under the accession numbers listed in Supplementary Table 1. Total genomic DNA (gDNA) extraction of the two specimens collected by us was performed by using the cetyl trimethyl ammonium bromide (CTAB) method [33]. To ensure that each DNA sample was of adequate quality for PCR, UV spectrophotometry (Nanodrop 2000, Thermo, Wilmington, DE) was used to determine the concentration and purity of DNA.
High-quality gDNA was then sent to the Novogene Co., Ltd. (Tianjin, China) for library preparation and genome sequencing. Short insert libraries (350 bp) were constructed using the NEB-Next Ò Ultra DNA Library Prep Kit for Illumina (NEB, USA). Whole-genome sequencing (WGS) was carried out on an Illumina Hiseq X Ten platform, generating 150 bp paired-end reads for each sample.

Mitogenome assembly and annotation
Approximately 5 Gb of raw data were obtained through WGS, and were then trimmed of adapters and low-quality reads (ratio of the bases less than Q20 > 30% or containing undetermined bases) by using the fastp v0.13.1 [34]. FastUniq v1.1 [35] and Musket v1.1 [36] were successively used for duplicate removing and error correction. Cleaned paired-end reads were de novo assembled using the SPAdes v3.14.1 software with k-mers of 21, 33, 55, 77, 99 and 127 [37]. Mitochondria-related contigs were identified and pooled for each sample by performing BLASTn search against the reference mitogenomes from B. cookei [38] and B. sorokiniana [39], which are closely related with the genera Exserohilum. We then filled gaps between these contigs by using the MITObim v1.9 to build closed-circular mtDNA macromolecules of the two Exserohilum samples [40]. In addition, both the assembled mitochondrial sequences were verified by NOVOPlasty [41].
The obtained complete mitogenomes of E. turcicum and E. rostratum were annotated according to the methods previously described by Wu et al. [42]. In short, we initially used Mfannot [43] and MITOS [44] to predict the PCGs, introns, unidentified open reading frames (uORFs), and non-coding RNA (ncRNA) genes among the mitogenome sequences, both based on genetic code 4. Then, the NCBI Open Reading Frame Finder (https://www.ncbi. nlm.nih.gov/orffinder) were used to modify or predict the PCGs and uORFs, which were further annotated by BLASTp searches against the NCBI non-redundant protein sequence database [45]. We verified the intron-exon borders of PCGs by using the exonerate v2.2 software [46] with the B. cookei mitogenome as reference [38]. Secondary structures of the tRNA genes were predicted by MITOS with default parameters [44], and were redrawn in Adobe Illustrator CS6. Graphical maps of the two Exserohilum mitogenomes were drawn using the online program OGDraw v1.2 [47].
In addition, for further comparative analyses, we downloaded another 11 complete mitogenomes of Pleosporales species (including B. cookei and B. sorokiniana that were previously used as references) from the NCBI GenBank database. Their detailed information (such as accession numbers) is provided in Supplementary Table 2.

Analyses of sequences and repetitive elements
Strand asymmetries of the 13 Pleosporales mitogenomes were assessed based on the following formulas: AT skew = [A -T] / [A + T], and GC skew = [G -C] / [G + C] [48]. Mitochondrial genome collinearity analysis of the two Exserohilum species and three closely related Bipolaris species was conducted by using Mauve v2.4.0 [49]. Based on genetic code 4, the sequence manipulation suite [50] was used to analyze the codon usage frequency and preference of the two Exserohilum mitogenomes. The DnaSP v6.10.01 software [51] was used to calculate synonymous (Ks) and nonsynonymous substitution rates (Ka) for 12 core PCGs (atp6, cob, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5, and nad6) in all 13 acquired Pleosporales mitogenomes. Based on the Kimura-2parameter (K2P) substitution model, MEGA v6.06 [52] was used to calculate the genetic distances between each pair of the 12 core PCGs. We performed BLASTn searches [53] of the two Exserohilum mitogenomes against themselves to detect if there were interspersed repeats or intragenomic duplications of large fragments, based on an E-value of < 1e-10. Tandem Repeats Finder [54] was used to detect tandem repeats (length > 10 bp) within the two Exserohilum mitogenomes. In addition, the contribution rates of different genetic compositions to the mitogenome expansion/contraction of two Exserohilum species were calculated according to the formulas:

Comparative mitogenomic and intron analyses
A comparative mitogenomic analysis was performed to evaluate the variations and conservations between different Pleosporales mitogenomes in genome size, base composition, GC content, gene number, intron number, gene arrangement, and gene content. According to the method described by Férandon et al. [55], introns within the cox1 genes of the 13 Pleosporales mitogenomes could be divided into different position classes (Pcls). Taking the cox1 gene of B. cookei as reference [38], cox1 genes from the other 12 Pleosporales species were aligned by Clustal W [56] to detect the insertion sites of introns. Pcls were named according to their insertion sites in the corresponding reference sequence. The same Pcls from different species usually have high sequence similarity and contain homologous intronic ORFs.

Mitochondrial phylogenomic analysis
In order to investigate the phylogenetic positions of the two Exserohilum species, we constructed a phylogenetic tree composed of 98 ascomycetes based on the concatenated mitochondrial gene set (nad1-6, nad4L, cox1-3, cob, rps3, atp6, atp8 and atp9). Both Taphrina deforman and T. wiesneri from Taphrinomycetes were appointed as the outgroup [57]. The software MEGA v6.06 [52] was first used to align individual mitochondrial genes, and then the SequenceMatrix v1.7.8 [58] was used to concatenated them into a combined mitochondrial gene set. Partition homogeneity test was used to detect potential phylogenetic conflicts between different mitochondrial genes. PartitionFinder v2.1.1 [59] was applied to determination of best-fit models of phylogeny and partitioning scheme of the gene set. We performed phylogenetic analysis using both Bayesian inference (BI) and maximum likelihood (ML) methods [60]. The ML analysis was implemented in the RAxML v8.2.12 with 1,000 bootstrap replications under the substitution model GTRGAMMA [61]. The BI analysis was performed using the MrBayes v3.2.6 [62]. Four simultaneous Markov chains were run starting from a random tree for 2,000,000 generations and trees were sampled every 100 generations. The first 25% of samples were discarded as burn-in, and the remaining trees were used to calculate values of Bayesian posterior probabilities (BPPs) in a 50% majority-rule consensus tree. Tree rendering was carried out using the FigTree v1.4.4 (https://tree.bio.ed.ac.uk/software/figtree/).

Data availability
The complete mitogenomes of E. turcicum and E. rostratum were deposited in the GenBank database under the accession numbers OK381869 and OK377062, respectively; and their raw sequencing data were deposited in the Sequence Read Archive (SRA) database under the accession numbers SRR16781348 and SRR16781347, respectively.

Features, genetic compositions and PCGs of two Exserohilum mitogenomes
Mitochondrial genomes of the two Exserohilum species, E. turcicum and E. rostratum, were circularly assembled, with the total sizes of 264,948 bp and 64,620 bp, respectively (Fig. 1A). Their GC contents were very close, with the values of 29.76% and 29.00%, respectively (Supplementary Table 2). For the E. turcicum mitogenome, both its AT and GC skews were positive, while in E. rostratum negative AT skew and positive GC skew were presented. Intronic regions occupied the largest proportion for both the two Exserohilum mitogenomes, reaching 62.25% and 31.53%, respectively (Fig. 1B). In the mitogenome of E. turcicum, protein-coding sequence was the second largest region accounting for 19.26%, whereas it ranked the third in E. rostratum, accounting for 27.64%. Intergenic spacers accounted for 14.25% and 29.63%, respectively, of the two mitochondrial genomes. Their ncRNA genes (including tRNAs and rRNAs), took up the smallest proportions, reaching only 4.23% and 11.20%, respectively.
Both the two newly sequenced Exserohilum mitogenomes had 12 typical core PCGs (atp6, cob, cox1-3, nad1-6, and nad4L) and a conserved core gene coding for the putative ribosomal protein S3 (rps3). Particularly, the cox1 genes in the two mitogenomes were seamlessly adjacent to the cox2, just resembling a fused gene. The mitogenome of E. turcicum contained two pseudogenes, which were respectively part of the two core PCGs, atp6 and nad2 (Supplementary Table 3). In addition, 34 free-standing ORFs were detected in the E. turcicum mitogenome, of which 21 were identified to encode homing endonucleases (HEs) and the others had unrecognized functions, whereas in E. rostratum, there were only six free-standing ORFs with unknown functions. A total of 70 introns harbouring 126 intronic ORFs were detected in the core PCGs and rRNAs of the E. turcicum mitogenomes, while in E. rostratum, only 17 introns with 17 intronic ORFs were found within the core PCGs. These intron-encoding ORFs mainly produced LAGLI-DADG and GIY-YIG HEs (Supplementary Table 3). In E. turcicum, most of the ORFs in introns encoded LAGLIDADG endonucleases, which were three times as many as GIY-YIG endonucleases. Analogously, the intronic ORFs in E. rostratum encoded more LAGLI-DADG endonucleases, two times as many as GIY-YIG endonucleases.

rRNA, tRNA, and codon analyses
Both the above Exserohilum mitogenomes had two rRNA genes, namely the small subunit ribosomal RNA (rns) and large subunit ribosomal RNA (rnl) (Supplementary Table 3). Nine and eight introns were respectively detected in the rnl and rns of the E. turcicum mitogenome, whereas no intron was found in both the two rRNA genes of E. rostratum. Therefore, E. turcicum possessed much longer mitochondrial rnl and rns genes (19,676 bp and 10,448 bp) than E. rostratum (3,408 bp and 1,682 bp).
A total of 32 and 29 tRNA genes were respectively identified in the mitogenomes of E. turcicum and E. rostratum, with lengths ranging from 70 bp (trnC-2) to 85 bp (trnY and trnS) (Fig. 2), which could encode 20 standard amino acids. Five (trnC, trnS, trnL, trnP and trnV) and four (trnL, trnP, trnN and trnS) tRNA genes were noticed to be duplicated in the E. turcicum and E. rostratum mito-genomes, respectively. Four copies of trnR and three copies of trnM were present in both the two mitogenomes. The copy number of trnN was three in the E. turcicum mitogenome. Overall, compared with the mitogenome of E. rostratum, the E. turcicum mitogenome contained three additional tRNA genes, i.e., trnN-3, trnC-2 and trnV-2. No introns had been found in these tRNA genes, and therefore all tRNA genes in the two mitogenomes could be folded into classical cloverleaf secondary structures. Of the 29 tRNAs shared by the two Exserohilum mitogenomes, only two variable sites were found in the trnP-1 and trnN-1, respectively.
ATG was the most commonly used start codon in the core PCGs of the two Exserohilum mitogenomes, except in the atp6 gene (TTG as the start codon) of E. turcicum. TAA was usually used as the stop codon in core PCGs of the two mitogenomes, followed by TAG. Codon usage analysis showed that the most frequently used codons in the two Exserohilum mitogenomes were TTA (for Leucine; Leu), AAA (for lysine; Lys), ATA (for isoleucine Ile), TCT (for serine; Ser), and AAT (for asparagine; Asn) ( Fig. 3 and Supplementary Table 4). The high frequent use of A and T in codons resulted in the high AT content of the two Exserohilum mitogenomes (average 70.62%).

Overlapping nucleotides and intergenic regions
There were six pairs of overlapping nucleotides detected in the mitogenome of E. turcicum, which were located across the neighboring genes orf139 and orf313 (À70 bp), orf70 and orf285 (À22 bp), orf130 and orf96 (À5 bp), orf329 and orf193 (À20 bp), orf555 and orf68 (À43 bp), and nad6 and orf388 (À85 bp), respectively. The E. rostratum mitogenome contained two pairs of overlapping nucleotides, of which one was located across the neighboring nad4 and orf99 genes (62 bp) and the other was located between orf99 and nad4L (1 bp) (Supplementary Table 3). A total of 37,766 bp and 19,147 bp of intergenic sequences were detected in the mitogenomes of E. turcicum and E. rostratum, respectively. Lengths of the intergenic spacers of the two Exserohilum mitogenomes ranged from 0 to 3,720 bp and from 0 to 1,446 bp, respectively. The longest intergenic segments in the E. turcicum and E. rostratum mitogenomes were located between orf168 and orf766 and between cox1 and orf138, respectively.

Analysis of repeat elements
Through BLASTn searches of the two Exserohilum mitogenomes against themselves, a total of 678 and 135 repetitive elements were separately detected in the mitogenomes of E. turcicum and E. rostratum (Supplementary Table 5). Sizes of these repeat regions ranged from 33 bp to 432 bp, with pairwise nucleotide similarities ranging from 81% to 100%. The longest repeat elements in the two mitogenomes were 432 bp and 424 bp, respectively. And they were observed to be located around the first exon of the cox1 gene and in the intergenic region between trnC and orf348, respectively. The overall length of the repetitive sequence in the E. turcicum mitogenome was 42,810 bp, accounting for 16.16% of its whole mitogenome, while the 7,012 bp of repetitive sequence in E. rostratum accounted for 10.83% of its mitogenome. In addition, a total of 50 and 26 tandem repeats were identified in the mitogenomes of E. turcicum and E. rostratum, respectively (Supplementary Table 6). Most of them in the two Exserohilum mitogenomes contained 2 to 4 copies, despite the highest copy number (65) observed in E. rostratum. The longest tandem repeat was found in E. turcicum, encompassing four repeat loci with a length of 93 bp. Tandem repeat sequences accounted for 1.43% and 3.00% of the mitogenomes of E. turcicum and E. rostratum, respectively.

Intron dynamics of the cox1 genes in Pleosporales
A total of 275 introns, most of them belonged to group Ⅰ, were detected in the mitogenomes of 13 Pleosporales species, each of which contained 1 to 70 introns (Supplementary Table 2). The great difference in intron number suggested that acquire/loss events of the introns had frequently appeared in the evolution of Pleosporales. In addition, the cox1 was estimated to be the largest host gene of introns in Pleosporales, containing 33.82% of the total introns, which was thus used for further analysis in the intron dynamic change among all the 13 Pleosporales species.
A total of 22 different Pcls were detected in the cox1 genes of 13 Pleosporales species (Fig. 4). The newly sequenced E. turcicum mtDNA was found to contain the largest number of Pcls (18) within the order Pleosporales, including two novel positions (P821 and P1307), whereas no introns were identified in the cox1 gene of Phaeosphaeria nodorum (Fig. 5). Half of the 22 Pcls were found to be widespread, among which P615, P731 and P1125, existing in nine of the 13 Pleosporales species, were the most common, followed by P1262 detected in seven Pleosporales species. P678 and P807 could be separately detected in only one of the 13 Pleosporales mitogenomes. Besides, the P678 of Stemphylium lycopersici seemed to be a species-specific Pcl that had been never detected in other ascomycetes.

Genetic distance and evolutionary rates of 12 core PCGs
As the rps3 gene was often absent from the Pleosporales mitogenomes, here only 12 core PCGs were used to calculate the genetic distances and substitution rates between each pair of the 13 Pleosporales species. It was found that the nad3 gene had the largest K2P genetic distance (average value 0.34), followed by the nad2 (average value 0.18), which meant that they exhibited the fastest mutation rates among the 12 core PCGs in Pleosporales (Fig. 6). The cox1 gene showed the lowest genetic differentiation between the 13 Pleosporales species, with an overall mean K2P distance of 0.10, indicating that this gene was highly conserved. Across the 12 core PCGs detected, the nad3 gene had the highest non-synonymous substitution rate (Ka) (average value 0.19) between the 13 Pleosporales species, followed by nad2 (average value 0.08), while the cox1 gene exhibited the lowest Ka value (average 0.02). The nad3 gene exhibited the highest synonymous substitution (Ks) rate (average value 3.73), while the lowest Ks value (average 0.78) among the 12 PCGs was observed in the cob gene. The overall Ka/Ks values for all detected core PCGs were<1, suggesting that these genes were subject to purifying selection.

Mitochondrial gene arrangement in Pleosporales species
Here, we compared the arrangements of mitochondrial core PCGs and rRNAs genes among the 13 Pleosporales species (Fig. 7), and found that the mitochondrial gene arrangement in Pleosporales often varied greatly, but was fairly conserved within genera, such as Exserohilum and Bipolaris. There were three Pleosporales mitogenomes that had no rps3 genes, including Phaeosphaeria nodorum, Didymella pinodes and Pithomyces chartarum. As mentioned earlier in the Exserohilum mitogenomes, an uninterrupted gene pair formed by the cox1 and cox2 was also found in all other Pleosporales species. The similar situation was observed between nad2 and nad3 and between nad4L and nad5. Besides, the three genes rnl, nad6 and rns together constituted a conserved gene block in most Pleosporales species except Shiraia bambusicola. Among the three closely related species of Bipolaris, no difference in gene order was found between B. sorokiniana and B. cookei, though position exchanges of two pairs of genes (cox1&cox2 and nad6&rns) were observed in the B. oryzae mitogen- ome. As for the two Exserohilum species, their mitochondrial gene arrangements were completely consistent, which appeared to be exactly the same as those of B. sorokiniana and B. cookei.

Mitogenome comparison
Comparative analysis showed that the mitogenome of E. turcicum was the largest one among all the 13 Pleosporales species available online (Supplementary Table 2). Its size (264,948 bp) was approximately 1.92 to 6.79 times larger than the other Pleosporales mitogenomes, indicating that E. turcicum might have experienced a huge mitogenome expansion during evolution. GC contents of the two Exserohilum mitogenomes (average 29.38%) were just slightly higher than those of other Pleosporales species (average 29.13%). Of the 13 Pleosporales species, each contained 15-49 PCGs (including core PCGs and free-standing ORFs) with the mean value of 20, which was close to the number of PCGs (19) in E. rostratum but much fewer than that (49) of E. turcicum. Besides, the E. turcicum mitogenome also had the largest numbers of introns (70) and intronic ORFs (1 2 6) among the 13 Pleosporales species, followed by the B. cookei mitogenome with 40 introns and 51 intronic ORFs. In addition, these 13 Pleosporales species gener-  Whole-mitogenome collinearity analysis was carried out between five closely related Pleosporales species, i.e., two Exserohilum and three Bipolaris fungi. A total of 11 homologous regions (A to K) were detected in each of the five Pleosporales mitogenomes (Fig. 8). Compared to the E. turcicum mitogenome, most of these homologous regions in the other four species were significantly reduced in size, especially for E. rostratum that seemed to have experienced a massive mitogenome contraction during evolution. Nevertheless, between the two Exserohilum mitogenomes, ten of the homologous regions (only except B) were arranged in the identical order. In addition, it was observed that, in the genus Bipolaris, the relative locations of most homologous regions were basically consistent as well. The homologous regions in B. sorokiniana and B. cookei had almost the same order (only except E), whereas the positions of five homologous regions (D, E, G, I, and J) varied between B. oryzae and the other two Bipolaris mitogenomes.
Remarkably, E. turcicum and B. cookei, belonging to two separate genera, had exactly the same order of homologous regions. In summary, the above Exserohilum and Bipolaris species, except B. oryzae, have a high degree of collinearity whether within or between genera.

Mitochondrial phylogenomic analysis
A phylogenetic tree involving 98 fungal species of Ascomycota was constructed based on both Bayesian inference (BI) and Maximum likelihood (ML) analyses of the 15 concatenated mitochondrial conserved PCG genes (Fig. 9). All major branches within this  tree were well supported (BPP ! 0.95; BS ! 89). Both T. deforman and T. wiesneri from Taphrinomycetes were appointed as the outgroup, while the other 96 species were well assigned into five clades in this topology, including Dothideomycetes, Sordariomycetes, Eurotiomycetes, Lecanoromycetes, and Pneumocystidomycetes (Supplementary Table 7). In the Dothideomycetes, different orders could be also well separated, namely Pleosporales, Cladosporiales, Mycosphaerellales, Dothideales, and Botryosphaeriales. In the order of Pleosporales, as expected, the two Exserohilum species tested here were clustered tightly with the highest support rates (BPP = 1; BS = 100%). And the genus Exserohilum was the most contiguous to the Bipolaris cluster containing three closely related species, B. cookie, B. sorokiniana and B. oryzae.

Discussion
Previous studies have indicated that fungal mitochondrial genomes were highly variable in size from 18.84 kb to 531.19 kb [24][25][26], and four main factors were reported to bring about this variation, including the introns, intergenic regions, repeat sequences, and plasmid derived dynamic change regions [63,64]. In this research, complete mtDNAs of two Exserohilum fungi were sequenced, assembled, and analyzed. Despite belonging to the same genus, these two species had vastly different mitogenome sizes, of which E. turcicum was the fourth largest (264,948 bp in length) among all reported mitogenomes of Ascomycetes, only smaller than the three macro-fungi Morchella crassipes (531,195 bp), Tuber calosporum (287,403 bp) and M. importuna (272,238 bp) [26,65,66]. Besides, it also became the largest one among all the 13 Pleosporales mitogenomes available online (Supplementary Table 1), which was about twice as big as the previously reported largest mitogenome of B. sorokiniana (137,775 bp) [39]. However, the size of the E. rostratum mitogenome (64,620 bp) was less than one-fourth of that of E. turcicum, even lower than the average value ($79 kb) of the other 11 Pleosporales species.
The earlier morphological identification, later multigene phylogenetic analyses based on nuclear loci, and current phylogenetic tree based on combined mitochondrial gene set (Fig. 9) all supported that the genera Exserohilum and Bipolaris had quite a close relationship [2,3]. Nevertheless, these two closely related taxa differed greatly in the mitogenome scale. Unlike Exserohilum, the three Bipolaris species used here, not surprisingly, possessed relatively similar mitogenome sizes, 124,887 bp, 135,790 bp and 137,775 bp, respectively, with a rangeability of only 10% ((maxmin)/min). Interestingly, mitogenomes of the two Exserohilum species were either much larger or far smaller than those of the genus Bipolaris, which implied that the two Exserohilum fungi might have undergone massive mitogenome expansion and contraction, respectively, during evolution. We found that, compared with E. rostratum and the three Bipolaris species, the mitogenome of E. turcicum had expanded by 310% (vs. E. rostratum), 112% (vs. B. oryzae), 95% (vs. B. cookei), and 92% (vs. B. sorokiniana), respectively. Of all the compositions of mitogenomes (described in the Results 3.1), intronic regions seemed to contribute the most to this expansion, with the contribution rates of 72.16%, 73.23%, 63.85% and 82.30%, respectively ( Supplementary Fig. 1), followed by the protein coding regions that contributed 16.56%, 14.67%, 27.78% and 18.87%, respectively. On the contrary, compared with E. turcicum and the three Bipolaris species, the mitogenome of E. rostratum had shrunk by 76% (vs. E. turcicum), 53% (vs. B. sorokiniana), 52% (vs. B. cookei), and 48% (vs. B. oryzae), respectively. Similarly, the intronic regions seemed to contribute the most to this contraction, with the contribution rates of 72.16%, 54.54%, 87.25% and 69.68%, respectively ( Supplementary Fig. 1). All these results indicated that the expansion/contraction of the Exserohilum mitogenomes was primarily due to their increase/decrease of intronic regions. Besides, repetitive sequences might also contribute to the expansion/contraction of the two Exserohilum mitogenomes, though their contribution rates were far below the intronic regions. In addition, plasmids and plasmid-like elements, which had been reported in some fungi as important factors leading to the change of mitogenome size [64,67,68], were not detected in the two Exserohilum mitogenomes.
Despite the huge distinction of the two Exserohilum mitogenomes in size, their core PCGs were very conserved in both number and sequence, which were consistent with the gene contents of other Pleosporales mitogenomes, consisting mainly of nad1-nad6, nad4L, cob, cox1-cox3 and atp6. However, there were another two core PCGs, atp8 and atp9, that could be commonly observed in most orders of Dothideomycetes, such as Mycosphaerellales [69], but were missing in all reported mitogenomes of Pleosporales including the genus Exserohilum [38,39,[70][71][72][73][74]. In addition, the two Exserohilum mitogenomes had a significant feature, i.e., the gene fusion between cox1 and cox2, which was also found in some other Pleosporales species, such as Corynespora cassiicola and S. lycopersici [73,74]. Estimation of genetic distances indicated that the sequences of cox1 genes were the most conserved between the 13 Pleosporales mitogenomes, whereas the nad3 genes exhibited considerable genetic differentiation. The Ka and Ks values varied not only between different core genes, but amongst the same genes of different species, especially for nad4 and nad4L, implying that in some Pleosporales species these genes may have evolved at a faster rate. All the 12 core PCGs of the 13 Pleosporales mitogenomes were identified to have been subject to purifying selection (Ka/Ks < 1), which meant that these genes had underwent a relatively conservative evolutionary process in functionality. In addition, the numbers of tRNA genes in the two Exserohilum mitogenomes (average 30.5) were just slightly higher than the average value (28.5) of other Pleosporales species. However, we noticed that, compared with E. turcicum, three tRNA genes (trnN-3, trnC-2 and trnV-2) were lost in the E. rostratum mitogenome; and besides, two variable sites occurred in the tRNAs trnP-1 and trnN-1 between the two Exserohilumum species, of which one was located in the acceptor stem (trnN-1) and the other was in the anticodon (trnP-1). The former may be a meaningless mutation, whereas the latter, if a non-synonymous mutation, can lead to the change of specific recognition of mRNA codon, and consequently affect the correct amino acid incorporation during translation [75,76]. The tRNA mutations that affect protein synthesis and metabolism have been found in some other eukaryotes [77,78].
Arrangement of the mitochondrial genes could provide a lot of reference information for understanding the phylogeny and evolution of eukaryotic species [79,80]. In animals, the mitochondrial gene rearrangement has been widely studied and a number of models associated with it have been proposed, such as the tandem duplication-random loss (TDLR) and duplication and nonrandom loss model [64,81]. By contrast, the rearrangement of fungal mitochondrial genes became more frequent, and their order and orientation could change in a larger scale [22]. In this research, a remarkable variation in gene arrangement was observed between different members of the Pleosporales, which was consistent with previous reports [72]. Nevertheless, several conserved gene blocks were found, like cox1-cox2, nad2-nad3, nad4L-nad5 and rnl-nad6-rns. Even between some genera, namely Exserohilum and Bipolaris, the gene arrangements were fairly conserved, which accorded with the fact that these two taxa were closely related and used to be classed as the same group. In addition, between B. oryzae and the other two Bipolaris species, obvious position exchanges of two pairs of mitochondrial genes (cox1&cox2 and nad6&rns) were observed, which led to a slight gene rearrangement within the genus. This limited exchange of gene positions between related species had been also found in the mitogenomes of Boletales [82]. Additionally, according to some previous studies [22,24], accumulation of the repetitive sequences in fungal mitogenomes correlated with the mitochondrial gene rearrangement. However, in this study, although a large number of repeat elements were detected in both the two Exserohilum mitogenomes, there was no gene rearrangement between them, which was consistent with the discovery in Boletus [82]. Overall, the mechanism of mitochondrial gene rearrangement in Pleosporales is quite important for understanding their evolution history, and deserves further investigation.
In eukaryotes, introns with various sizes and frequencies are present in nuclear, mitochondrial, and plastidial genes from different kingdoms [83]. Mitochondrial introns are classified in two main groups, namely group I and II, of which the former was mainly found in fungi and the latter was more common in plants [22]. In this study, a total of 275 introns, most of them belonged to group I, were detected in 13 Pleosporales mitogenomes, and each species had 21 introns on average. As more than a third of these introns were found in the cox1, this gene was thus selected to analyze the dynamic changes of mitochondrial introns among the 13 Pleosporales species. By looking into the 22 Pcls detected in the cox1 genes, we noticed that the 13 Pleosporales species had a great difference in the number and type of mitochondrial introns, which indicated that the gain/loss events of introns might have occurred in the evolution of Pleosporales. P615, P731 and P1125 were found to be the most conserved intron insertion sites in the order Pleosporales, and were therefore considered to be derived from the common ancestor of Pleosporales species. Per contra, each of the Pcls, P678, P807, P821 and P1307, only existed in one of the 13 Pleosporales mitogenomes. And these scarce introns had been detected in some distantly related species, such as Arthrobotrys musiformis [84] and Monilinia fructicola [85], suggesting potential horizontal gene transfer events. Although E. turcicum possessed the second largest number of introns (18) in the cox1 gene, only one less than that of Agaricus bisporus [55], there were no novel Pcls detected in it, whereas in S. lycopersici, a novel Pcl (P678) was probably found, which had been never detected in all other known Ascomycota mitogenomes. For fungal mitogenomes, group I introns were considered to be an important genetic element [82,83], and their dynamic changes could significantly affect the size and organization of mitogenomes [29,65].
Dothideomycetes is the largest class of Ascomycota, including a great number of economically important plant pathogens [86]. Apparent characteristics of many closely related species or species complex in the Dothideomycetes are easy to be confused, which makes it difficult to distinguish these taxa accurately only according to the morphology. Therefore, it is very urgent to find reliable molecular markers to solve the above classification problem. Although some researchers had begun to use single-gene and multi-gene phylogenetic analyses to explore the genetic relationships within the class Dothideomycetes [87,88], it is still often a hard task to identify some complicated genera and species. In recent years, due to their advantages over nuclear genomes including the uniparental inheritance and accelerated evolutionary rate [89], mitogenomes have been widely used in phylogeny and population study of plants, animals, and some kinds of fungi [22,24,64,82]. Here, we established a highly supported phyloge-netic tree of 98 Ascomycota species based on BI and ML analyses of the combined mitochondrial gene set. All these selected species were well gathered or divided into different independent clades and subclades, each of which exactly corresponded to one fungal taxa (such as class and order). The two Exserohilum species were observed to be the most closely related to three Bipolaris species, which was consistent with the previous analyses of multigene phylogeny based on nuclear loci [3]. Our result indicated that mitochondrial genomes might be useful markers (such as the cox1, cox2, and cob genes) for fungal phylogenetic analysis in different taxon level. These mitochondrial barcoding genes have been also used to study biological evolution of animals [90,91]. However, the number of known mitogenomes in Ascomycota is still limited, which restricts our more accurate reconstruction of the phylogenetic relationships between related species. A better understanding of mitogenome evolution and phylogeny within Ascomycota would require the sequencing of more members from this group.

Conclusion
In this study, the mitochondrial genomes of two Exserohilum species, E. turcicum and E. rostratum, collected from China were sequenced, assembled, annotated and compared. The mitogenome of E. turcicum was the largest (264,948 bp) among all reported Pleosporales mitogenomes, whereas the E. rostratum mitogenome was extraordinarily small (64,620 bp), actually less than a quarter of the former. Through the comparative analysis of genetic compositions, increase/decrease of the intronic regions was considered as the main contributing factor of their mitogenome expansion/contraction. Despite the huge difference in size between the two Exserohilum mitogenomes, we confirmed that they shared the same set of 13 core PCGs and two rRNAs in identical order and orientation. All the 12 core PCGs (excluding rps3) of 13 Pleosporales mitogenomes were identified to have been subject to purifying selection (Ka/Ks < 1). In addition, a total of 93 group I introns, classified into 22 Pcls, were differentially distributed in the cox1 genes of 13 Pleosporales species, of which E. turcicum possessed the second largest number (18). Loss/gain and horizontal transfer events of these introns might promote size variation of the Pleosporales mitogenomes. Mitogenome-wide phylogeny well exhibited the genetic relationships between selected ascomycetes including the two Exserohilum species, suggesting that mitochondrial genes were reliable molecular markers for phylogenetic analysis of Ascomycota. This study served as the first investigation on Exserohilum mitogenomes, which will lay the foundation for further understanding the genetic evolution, species differentiation, and ecological adaptation of Exserohilum species.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


The immunology of low-grade gliomas

Low-grade gliomas (LGGs), which harbor an isocitrate dehydrogenase (IDH) mutation, have a better prognosis than their high-grade counterparts; nonetheless, they remain incurable and impart significant negative impacts on patients’ quality of life. Although immunotherapies represent a novel avenue of treatment for patients with LGGs, they have not yet been successful. Accurately selecting and evaluating immunotherapies requires a detailed understanding of LGG tumor immunology and the underlying tumor immune phenotype. A growing body of literature suggests that LGGs significantly differ in their immunology from high-grade gliomas, highlighting the importance of investigation into LGG immunology specifically. In this review, the authors aimed to discuss relevant research surrounding the LGG tumor immune micro environment, including immune cell infiltration, tumor immunogenicity, checkpoint molecule expression, the impact of an IDH mutation, and implications for immunotherapies, while also briefly touching on current immunotherapy trials and future directions for LGG immunology research.

L ow-grade gliomas (LGGs) are slow-growing glial tumors with the potential for malignant transformation. Their incidence is highest in patients in their 3rd and 4th decades of life, and both the tumor and the treatments can significantly affect the quality of life of patients. Seizure is the most common presenting symptom in the majority of patients, although headache and focal neurological deficits are also frequently seen. 1,2 Indeed, LGGs are more often present in highly functional areas of the brain, with the insula and supplemental motor area being the two most commonly afflicted locations. 3,4 The standard of care for LGGs includes surgical intervention, with the goal of providing maximal safe resection and acquisition of adequate tissue for detailed molecular and genetic characterization. 5 This is becoming increasingly important, as the molecular subtypes of LGGs and their unique characteristics are becoming better defined. 6 The utility of molecular and genetic characterization for LGGs is most evident in the latest 2021 WHO guidelines, which place an increased emphasis on molecular diagnostics. 6 These updated guidelines for gliomas focus specifically on mutations in the isocitrate dehydrogenase 1 (IDH1) and IDH2 genes, which have been associated with significantly improved outcomes in patients. 7,8 As a result, adult LGG is now, by definition, IDH mutant, including both astrocytoma and oligodendroglioma (which also harbor a 1p/19q-codeletion). 6 This is partly due to the differences in the natural history and outcomes of patients with IDH-wild-type and IDH-mutant gliomas. 8 Few effective new therapeutic strategies for gliomas have been developed in the past decade. Immunotherapy, which targets the immune system against tumor cells, has been successful in other cancer types and is a major area of research for gliomas, albeit predominantly for highgrade gliomas (HGGs). Immunotherapies that are under investigation for gliomas include chimeric antigen receptor (CAR) T cells, systemically delivered checkpoint inhibitors (CPIs), oncolytic viral therapies, other viral gene therapies, and cancer vaccines, all of which have unique limitations. 9 While previously thought to be immune-privileged, the brain has been highlighted in recent studies as being accessible to the immune system, suggesting that immunotherapies may hold promise in treating CNS tumors, including glioblastoma. 10 Although initial trials of immunotherapies in HGGs have been met with limited success, 11,12 literature supporting different immunology in LGGs, and the hypothesis that altering the immune environment is a gradual process that might be more effective for LGGs than HGGs, have led to excitement about immunotherapies for treatment of LGGs. In this review, we summarize the immune phenotypes of patients with LGGs and how this may impact immunotherapies for LGGs as well as challenges and future directions in LGG immunology.

Immune Cell Infiltration
Immune cell infiltration of tumors has been shown to play a role in tumor progression and has been used for cancer prognostication. [13][14][15] Immune cell infiltration can also provide insight into how to best orchestrate a robust antitumor immune response using immunotherapy. Assessment of immune cell infiltration can be done using a variety of techniques including bulk and single-cell RNA sequencing, bulk genomic DNA samples, flow cytometry, cytometry by time of flight, immunofluorescence, and additional methodologies.
In a study of microarray and bulk RNA sequencing data from 2249 gliomas, Zhang et al. demonstrated that reduced tumor purity (i.e., increased invasion of a glioma by nontumor cells such as immune and stromal cells) was correlated with increased malignancy and reduced survival. 16 The authors also showed that tumors with reduced glioma purity had increased immune cell scores, especially macrophage and neutrophil populations, relative to higher purity samples. Similarly, Aran et al., in an analysis of multiple cancer types, demonstrated an inverse relationship between tumor purity and mutational burden, with LGGs as one of the lowest mutational burdens and highest tumor purities in their analysis. 17 These findings on tumor purity suggest that LGGs are relatively immune and stromal cell sparse relative to HGGs and other cancer types.

T-Cell Compartment
More detailed analyses of the LGG immune phenotype have shed additional light on specific immune cell infiltration. A study by Berghoff et al., which focused on the T-cell compartment in IDH-mutant and IDH-wild-type gliomas using immunohistochemistry, found that tumorinfiltrating lymphocytes and programmed death-ligand 1 (PD-L1) expression were significantly lower in IDHmutant gliomas compared with IDH-wild-type tumors. 18 The authors also found that LGGs had significantly higher levels of PD-L1 methylation, likely explaining the reduced expression. Additional studies by Klemm et al. and Friebel et al. leveraged multifaceted approaches to assessing the CNS tumor microenvironment and similarly demonstrated low tumor-infiltrating lymphocyte levels in LGG relative to glioblastoma and brain metastasis. 19,20 These findings suggest that LGGs, HGGs, and brain metastases lie on a spectrum of immunogenicity and T-cell infiltration, with brain metastases demonstrating the most T-cell-inflamed immune phenotype and LGGs the least-inflamed phenotype, with HGGs falling in between the two. Interestingly, similar trends have also been seen in the context of regu-latory T cells (Tregs), whereby Treg infiltration positively correlated with the glioma tumor grade. 20,21 Friebel et al. also showed higher Treg infiltration in brain metastasis compared with HGGs and LGGs, similar to the trend seen with T cells in general. The differences in T-cell infiltration between LGGs, HGGs, and brain metastases suggest a tumor-intrinsic process driving the surrounding immune phenotype, rather than just the location intracranially.
The underlying biology of the differences in the T-cell compartment between various CNS tumors, especially IDH-mutant and IDH-wild-type gliomas, is an area of active research. Chongsathidkiet et al. implicated tumorimposed loss of S1P1 from the T-cell surface in sequestration of T cells in the bone marrow in multiple intracranial murine tumor models, including both glioma and metastasis, suggesting a common immunosuppressive pathway for brain tumors. 22 Kohanbash et al., however, sought to investigate the specific role of IDH mutations in shaping the tumor microenvironment to further understand the unique immunophenotypes in IDH-mutant and IDHwild-type gliomas. 23 Using The Cancer Genome Atlas (TCGA), the authors found reduced expression of T-cell lymphocyte-associated genes as well as CXCL10, which is implicated in the trafficking of CD8 T cells into brain tumor sites, in IDH-mutant tumors compared with IDHwild-type tumors. They then introduced mutant IDH1 or the downstream metabolite of mutant IDH, oncometabolite R-2-hydroxyglutarate (2HG), to syngeneic mouse glioma models. Expression of 2HG or expression of mutant IDH1 reduced CXCL10 levels secondary to decreased production of STAT1 (a CXCL10 regulator); this was also associated with reduced T-cell infiltration in implanted tumors. Finally, IDH-C35, an inhibitor of mutant IDH1, reversed the immunophenotype seen in mutant tumors and increased the efficacy of vaccine immunotherapy. 23 These findings implicated the production of 2HG by mutant IDH1 in molding the eventual immunophenotype seen in these tumors. Interestingly, mutant IDH does not seem to have as strong of an effect on CD4 T-cell levels, although previous research has demonstrated that CD4 T cells are, in fact, capable of generating an immune response against the R123H IDH mutation. 24,25

Myeloid Cell Compartment
While T cells play an undeniable role in generating a robust antitumor immune response, significant consideration should also be paid to the myeloid compartment of a tumor's immune phenotype. Myeloid cells, which make up a much larger percentage of the immune cell compartment relative to T cells, can play a range of roles in supporting or hampering an antitumor immune response, and differences in the infiltration of myeloid cells can have an important impact on the efficacy of immunotherapies ( Fig. 1). Key classes of myeloid cells in the tumor microenvironment include dendritic cells, neutrophils, tumor-associated macrophages (TAMs; including monocyte-derived macrophages [MDMs] and microglia), and myeloid-derived suppressor cells (MDSCs). [26][27][28] Indeed, TAMs comprise the majority of IDH-mutant and IDHwild-type tumors, significantly outnumbering the number of other immune cell types, including T cells. 19,20 Interest-ingly, the dominant type of TAM differs between IDHmutant and wild-type tumors, with IDH-mutant TAMs being composed primarily of microglia and IDH-wild-type tumors (and brain metastases) associated with primarily MDMs. 19,20,29 The microglia present in IDH-wild-type tumors and brain metastases also appear to be more of an activated phenotype than those found in IDH-mutant tumors, as determined by differences in CD14 and CD64 genes and cell morphology. 20 Within LGGs, Müller et al. demonstrated significantly reduced infiltrating microglia in oligodendrogliomas relative to astrocytomas; LGGs with higher MDM infiltration also had significantly reduced survival. 29 The difference in microglial and MDM cells in IDH-mutant astrocytomas and oligodendrogliomas was also noted by Venteicher et al. 30 In their study, astrocytomas also demonstrated increased MDMs as the tumor grade increased, which corresponded with an increase in the expression of endothelial-specific genes. The association between MDMs and endothelial genes may be due to increased access of MDMs to the tumor microenvironment from systemic circulation. It is also important to note that tumor-infiltrating MDMs lie on a spectrum; Venteicher et al. noted that MDMs in the glioma microenvironment expressed less macrophage-specific genes than those found in brain metastasis, suggesting a range of macrophage states rather than clearly defined and circumscribed cellular states. 30 Given the difficulty in distinguishing MDSCs from other myeloid subtypes in the tumor microenvironment, there is a paucity of data on MDSC infiltration in LGGs; they are likely being labeled as either MDMs or neutrophils (depending on if they are granulocytic or monocytic) in studies to date.
The intratumoral forces driving the patterns seen in myeloid cell infiltration in LGGs are still unclear. As previously discussed, the association between tumor grade, endothelial gene expression, and MDM infiltration noted by Venteicher et al. may suggest that increased exposure to the systemic circulation, through increased tumor vas- cularity and blood-brain barrier breakdown, may contribute to differences in myeloid infiltration between LGGs and HGGs. 30 Increased MDMs in low-grade astrocytomas relative to oligodendrogliomas also potentially indicate a different tumor cell-intrinsic mechanism. Seeking to evaluate the effect of tumor genotype on the myeloid cell compartment, Friedrich et al. assessed alterations in the immune phenotype of tumors in response to IDH mutation using a variety of molecular tools. 31 Similar to the aforementioned studies, the authors determined that microglia and macrophages demonstrated wide separation on principal component analysis corresponding with IDH status following RNA sequencing in HGGs. Using the GL261 mouse model of glioblastoma, they found that the myeloid compartment of the tumor microenvironment was significantly impacted by tumor IDH status. Mice implanted with IDH-mutant tumors demonstrated lower expression of the CCL2 chemokine (which recruits myeloid cells), less immunoactivated microglia, and delayed MDM infiltration relative to wild-type tumors. Interestingly, MDMs from IDH-mutant tumors showed an increased ability to suppress T-cell production of interferon-γ (IFNγ) and incite expression of the programmed cell death protein 1 (PD-1) exhaustion marker. Inhibiting IDH led to a proinflammatory phenotype in macrophages, while the exogenously administered 2HG recapitulated the MDM (but not microglia) findings seen in IDH-mutant tumors, directly implicating 2HG in the immunosuppressive macrophage phenotypes seen in IDH-mutant tumors. The authors then further identified tryptophan metabolism, specifically the enzyme tryptophan 2,3-dioxygenase (TDO) and its product kynurenine, as playing a key role in aryl hydrocarbon receptor (AhR) activation and the downstream immunosuppressive effects previously characterized. 31 These results provide some insight into how IDH-mutant tumors, which, by definition, are LGGs, suppress the surrounding immune microenvironment.
A recent study by Ludwig et al. identified an alternative pathway in which IDH-mutant gliomas can generate an immunosuppressive phenotype both locally and systemically through tumor-derived glioma small extracellular vesicles (TEX). 32 The authors found that IDH-mutant gliomas produced larger, more immunosuppressive, and morphologically distinct TEX from IDH-wild-type tumors. Injection of IDH-mutant TEX into wild-type tumors led to reductions in tumor-infiltrating natural killer cells and increases in immunosuppressive Tregs and mono-MDSCs. Interestingly, injection of IDH-wild-type TEX into IDH-mutant tumors increased effector lymphocytes and M1 macrophages and reduced MDSCs. 32 The studies by Ludwig et al. and Friedrich et al. highlight the still-developing theories on how IDH impacts glioma tumor immunology. 31,32 Additional studies are required to fully understand the impact of IDH on immune cell infiltration, specifically in LGGs, and potential therapeutic vulnerabilities.

Antigen Presence and Presentation
Tumor mutational burden is seen as a potential indicator of immunogenicity, as in many cancer types it is thought that an increased mutational burden is also likely indicative of the overall tumor antigen number, providing more "targets" for the immune system. [33][34][35][36][37] This was demonstrated in a study by Samstein et al., who assessed the tumor mutational burden of 1662 patients with CPI treatment and 5371 without CPI. 37 The authors found that, for most cancers, patients with a higher mutational burden had better responses to CPI treatment, despite inconsistent mutational cutoffs between cancer types. Interestingly, gliomas did not follow the same trend as other cancers, trending instead toward worse survival in those patients with high tumor mutational burdens. Alternatively, hypermutated gliomas have been shown to have low T-cell infiltration, high intratumoral heterogeneity, and a poor response to anti-PD-1 therapy. 38 This may be due to temozolomide-induced mutations, which are frequently single-nucleotide variant mutations and do not adequately induce an antitumor immune response. In fact, in HGGs, low mutational burden has been associated with improved outcomes following viral or CPI therapy. HGGs with a lower mutational burden had increased inflammatory gene signals, suggesting that neoantigens were removed from the tumor via immunoediting. 39 In LGGs, a higher mutational burden has similarly been associated with worse outcomes; 40 this is especially true if patients develop temozolomide-related hypermutation upon transformation to a higher grade. 41 How hypermutation specifically impacts the immune tumor microenvironment in LGGs is still not fully elucidated and warrants additional investigation.
Similarly, there is a paucity of literature on antigen presentation and dendritic cells in LGGs. Dendritic cells are critical for immune surveillance of the CNS. They can capture antigens, process them, and present them to T cells (both cytotoxic and helper) and B cells, thereby acting as a bridge between the innate and adaptive immune systems. The literature on dendritic cell infiltration in LGGs is limited, but there appear to be few dendritic cells in both HGG IDH-wild-type and LGG IDH-mutant tumors. 19 Interestingly, it appears that LGG IDH-mutant tumors have a reduced capacity for antigen presentation, as IDH-wildtype tumors have increased infiltration of MDMs expressing major histocompatibility complex (MHC) class II human leukocyte antigen-DR isotype (HLA-DR). HGG IDH-wild-type MDMs also have increased expression of antigen presentation gene sets, antigen processing associated pathways, and MHC class I presentation gene sets. 19 It must be noted, however, that HGG IDH-wild-type tumors remain resistant to immunotherapy and MDMs have been implicated in tumor progression, making it difficult to interpret their possibly increased ability for antigen presentation. Additional analysis is needed to compare the functional antigen presentation ability of MDMs from IDH-mutant tumors with those from IDH-wild-type tumors to better understand their ability and how to harness them for immunotherapy.

Checkpoint Molecule Expression
Checkpoint molecules are frequently pathologically coopted by cancers to avoid immune surveillance and killing. This has driven the development of monoclonal antibodies to block checkpoint molecules, also known as a CPI. While checkpoint molecule expression in the tumor microenvironment may not always correlate with tumor responsiveness to checkpoint inhibition, it is still important to consider. Garber et al. sought to assess PD-1 and PD-L1 expression in various glioma grades. 42 The authors found that both PD-1-positive tumor-infiltrating lymphocytes and PD-L1 expression in tumor cells were correlated with increasing tumor grade. 42 Interestingly, there did not appear to be a relationship between IDH-mutant tumors and either PD-1-positive tumor-infiltrating lymphocytes or PD-L1-expressing tumor cells. 42 A subsequent study by Berghoff et al. found significantly higher levels of PD-1-positive tumor-infiltrating lymphocytes and PD-L1 expression in IDH-wild-type gliomas compared with IDH-mutant gliomas. 43 They also found that IDH-mutant tumor cells had higher PD-L1 promoter methylation, perhaps explaining their observations. The findings by Berghoff et al. correlate with a study of 1024 glioma samples by Liu et al. that demonstrated higher CTLA-4 expression in higher-grade IDH-wild-type tumors compared with lower-grade IDH-mutant tumors. 44 Overall, it appears that IDH-wild-type HGGs express higher levels of checkpoint molecules, potentially making them more attractive candidates for CPI immunotherapy. However, further investigation is needed to elucidate the underlying molecular biology informing the observed trends.

Current Immunotherapy Trials for LGGs
As discussed in this review article, LGGs have unique tumor immune phenotypes, with their IDH mutations potentially making them more or less vulnerable to specific types of immunotherapies. Compared with HGGs, there is a relative paucity of clinical trials evaluating immunotherapies for IDH-mutant LGGs. This is likely due to a variety of factors, including the better prognosis associated with LGGs and their lower incidence, resulting in the need for relatively longer clinical trials to evaluate the efficacy of a novel treatment. An additional difficulty in developing immunotherapies for LGGs is the differences in immune gene expression between HGGs and LGGs, which can hinder the ability for HGG immunotherapy trials to be translated to LGGs. For example, CAR T-cell targets for HGGs, such as EGFRvIII and IL13Ra2, are not found in LGGs, preventing their use in LGGs and requiring the identification of LGG-specific targets.
However, immunotherapeutic modalities, such as vaccines, may offer a safe and effective option for patients with LGGs due to the slower growth rate of the tumors (in contrast with HGGs), which should allow sufficient time for multiple immunizations and higher levels of antiglioma immunity. Furthermore, the generally mild toxicity of vaccines may have advantages over chemotherapy or radiation therapy for long-term cognitive and quality-of-life impairments. Okada et al. conducted the first vaccine trial that was specifically designed for patients with LGGs and reported that these patients, especially those who had not received chemotherapy or radiation therapy, demonstrated an excellent immunological response to the vaccine. 45 Immunotherapies that are being evaluated in IDH-mutant LGGs now include CPIs; polyinosinic-polycytidylic acid stabilized by lysine, and carboxymethylcellulose (poly-ICLC); cancer vaccines; and IDH-mutant inhibitors. A list of immunotherapy clinical trials including IDH-mutant LGGs and their underlying rationales can be found in Tables 1 and 2, respectively. Thus far, no immunotherapy has demonstrated a robust survival impact in clinical trials for LGGs, although the slow rate of progression in these gliomas means that very few trials have completed their evaluation (Table 1). Promising newer developments in immunotherapy for LGGs include the combination of IDH-mutant inhibitors with another immunotherapy such as a CPI. This may help to reduce the immune suppression associated with IDH-mutant tumors before adding an additional agent to stimulate the immune system. 46

Conclusions
Despite their slow-growing nature, LGGs lead to significant morbidity and mortality. They are, in part, defined by the presence of mutant IDH, which may alter their underlying biology in a variety of ways compared with IDH-wild-type tumors. Recent studies using novel multifaceted immune-profiling approaches have significantly increased our understanding of LGGs, which have tumor immune microenvironments that are uniquely immunosuppressed compared with IDH-wild-type HGGs and brain metastases; LGGs have been shown to have reduced T-cell infiltration and decreased infiltration of myeloid cells with the potential for antigen presentation. Initial studies have demonstrated the importance of mutant IDH and its byproduct 2HG in mediating this unique immune phenotype, but a complete understanding of the molecular pathways implicated is lacking. As a result, significant investigation is needed to better understand the LGG immune microenvironment, underlying biological pathways, and implications for immunotherapies.
There are a number of areas within LGG research that warrant additional investigation and development. One such area of improvement that will allow for future advancement in the realm of immunology and immunotherapy for patients with LGGs is the generation of preclinical cell lines that are better biomimetic models for human LGG tumors. Currently, the majority of mouse glioma models are created in an effort to recapitulate glioblastoma rather than LGGs. 47 Indeed, studies investigating how IDH mutation affects the immune microenvironment have primarily been performed in murine glioblastoma models. Introducing IDH mutations into murine glioblastoma cell lines likely does not accurately model the tumorigenesis and pathology of human LGGs, given the significant differences in tumor invasion and breakdown of the bloodbrain barrier, which may impact immune profiles. An additional difficulty with investigating immunotherapy for LGGs is the longer and variable patient survival associated with these tumors, which makes clinical trials logistically challenging and more expensive compared with glioblastoma trials. Finally, the systemic impact of gliomas on the immune system is better defined for HGGs than LGGs. 22 Although our previous vaccine trial in patients with LGGs suggested that those patients can mount superior immune responses against peptide vaccines compared with patients with HGGs or pediatric patients with gliomas, ad-ditional insight into how LGGs impact systemic immunity is needed. 45 Whether or not LGGs lead to the T-cell sequestration similar to that seen in HGGs is unclear. This is especially relevant because patients with LGGs may be treated with temozolomide for long periods of time, which can also lead to systemic immune suppression and might interfere with certain immunotherapy treatments. 48  Seek to train & target the immune system against tumor-associated antigens to stimulate an antitumor immune response IDH-mutant inhibitors Variety of mechanisms; may contribute to increased antitumor immune activity by reversing or limiting IDH-mutant-associated immune suppression