Systematic analyses with genomic and metabolomic insights reveal a new species, Ophiocordyceps indica sp. nov. from treeline area of Indian Western Himalayan region

Ophiocordyceps is a species-rich genus in the order Hypocreales (Sordariomycetes, Ascomycota) depicting a fascinating relationship between microbes and insects. In the present study, a new species, Ophiocordyceps indica sp. nov., is discovered infecting lepidopteran larvae from tree line locations (2,202–2,653 m AMSL) of the Kullu District, Himachal Pradesh, Indian Western Himalayan region, using combinations of morphological and molecular phylogenetic analyses. A phylogeny for Ophiocordyceps based on a combined multigene (nrSSU, nrLSU, tef-1α, and RPB1) dataset is provided, and its taxonomic status within Ophiocordycipitaceae is briefly discussed. Its genome size (~59 Mb) revealed 94% genetic similarity with O. sinensis; however, it differs from other extant Ophiocordyceps species based on morphological characteristics, molecular phylogenetic relationships, and genetic distance. O. indica is identified as the second homothallic species in the family Ophiocordycipitaceae, after O. sinensis. The presence of targeted marker components, viz. nucleosides (2,303.25 μg/g), amino acids (6.15%), mannitol (10.13%), and biological activity data, suggests it to be a new potential source of nutraceutical importance. Data generated around this economically important species will expand our understanding regarding the diversity of Ophiocordyceps-like taxa from new locations, thus providing new research avenues.


. Introduction
Mushrooms are widely accepted as macrofungi with distinctive, immense, easily notable, and collectable fruiting bodies, serving distinct ecological and economic functions. However, few mushrooms have been reported as toxic (Horowitz, 2015;Govorushko et al., 2019;Tran and Juergens, 2022) due to the presence of harmful secondary metabolites and heavy metals, especially arsenic (Li Y. et al., 2019). Many new species of medicinal mushrooms have been identified, based on morphological characteristics combined with phylogenetic approaches. Ancient oriental customs have emphasized the significance of different medicinal mushrooms with prominent bioactive components which are being valued throughout the world. Among these, the genus Ophiocordyceps is one of the most fascinating insectfungal associations deciphering immense nutraceutical and therapeutic potential.
To date, among all the studied Ophiocordyceps spp., O. sinensis (Berk.), synonym Cordyceps sinensis, commonly known as "Himalayan Viagra" (Yarsagumba, Winter Worm Summer Grass) and Dong Chóng Xià Cǎo in Chinese, is world's most expensive and extremely rare natural medicinal resource. It has a very restricted habitat in the harsh environment of alpine meadows (3,000-5,500 m above mean sea level, AMSL) of China, Tibet, Bhutan, Nepal, and India (Li et al., 2011). As a result, harvesting this rare and valuable mushroom is highly challenging and even dangerous as the collectors have to travel through tough terrains. In India, it is known as Keeda Jari or Keeda Ghas and is reported from Uttarakhand (U.K.), Sikkim, and Arunachal Pradesh . Ophiocordyceps sinensis is commonly recognized as "soft gold" (Woodhouse et al., 2014), "fungal gold" (Panicker, 2017), or "organic gold" (Pouliot et al., 2018), and reported to be even pricier than gold (Li et al., 2015). Owing to the findings of bioactive components in relation to anti-aging, anti-tumor, antioxidant, antiinflammatory, and immuno-modulatory activities, Ophiocordyceps has attracted substantial interest in the modern medicine sector, in addition to its traditional consumption (Lou et al., 2019;Li Y. et al., 2021;Zhang et al., 2021). There was a 900% rise in its price between 1998 and 2008, and top-quality Ophiocordyceps were sold at $145,000/kg Li Y. et al., 2021). In addition, its high price and demand led to overharvesting, overexploitation, damage to its fragile mountain environment, and a sharp drop in its yield, which brought it close to extinction (Hopping et al., 2018). Therefore, conservation, sustainable harvesting, and identifying its alternatives are strongly advocated and attempts were made for its in vivo cultivation (Zhou et al., 2014;Qin et al., 2018;. Nevertheless, in vitro cultivation of its anamorph viz. Hirsutella sinensis and its associated fungus, Samsoniella hepiali (syn. Paecilomyces hepiali), have been achieved by various researchers. As the research into Ophiocordyceps has progressed, new products, such as capsules, powder, and solutions, have been introduced in the global market (Zhou et al., 2009;Au et al., 2012;Hopping et al., 2018;Yang et al., 2018). Some of the products are Pure Cordyceps TM (Aloha Medicinals Inc.), Jin Shui Bao (Jiangxi Jinshuibao Pharmaceutical Co. Ltd.), Cordyceps capsule (Geneferm Biotechnology Co. Ltd.), and Bailing capsules (Hangzhou East China Pharmaceutical Co., Ltd.; Zhou et al., 2009;Elkhateeb and Daba, 2022). In India, wild O. sinensis and some mycelial-based products are being used as nutraceuticals and immunity boosters (Sharma, 2004;Panda and Swain, 2011;Maity, 2013). However, the demand, market, and price for these products are cheaper than the wild intact Ophiocordyceps.
Keeping in view the challenges and potential of O. sinensis, the identification, and positioning of any Ophiocordyceps species from the wild inhabiting low-lying and approachable areas is the need of the hour. Under these circumstances, the present study was planned, and exploratory surveys were conducted between 2015 and 2021 in different districts (Kangra, Chamba, Mandi, Kullu, and Lahaul and Spiti; 1,000-4,800 m AMSL) of Himachal Pradesh (H.P.) India, to explore the presence of any Ophiocordyceps species of commercial importance. During our surveys, one kind of unknown pathogenic fungus infecting larvae of Thitarodes sp. buried in soil was collected from tree line locations of Kullu, H.P., Indian Western Himalayas. The present study was conducted with the objective to confirm the taxonomic position of the collected specimens, and we introduced a new taxon of Ophiocordyceps of commercial importance from low-height areas (2,202-2,653 m AMSL) of Himachal Pradesh, India, with morphological descriptions coupled with multigene phylogeny, genomic, and metabolomic estimations.    indica). Samples were observed and collected for the first time at 2,479 m AMSL (Figures 2A-D) in July 2016. In subsequent years, specimens were further recorded and studied in almost similar habitats at different locations and altitudes (2,202-2,653 m AMSL) in the same district ( Figure 3A; Supplementary Videos 2, 3). Interestingly, the morphology and habits of these specimens were somewhat similar to O. sinensis

. . Morphological studies
The specimens were collected along with their host, the larva of Lepidoptera, and were kept in plastic bags and air-tight containers. Macromorphological characteristics were recorded during the collection time. The specimens were preserved by air drying followed by deep freezing and were stored in the Entomology Laboratory, CSIR-IHBT with voucher numbers. Furthermore, the specimens were studied and identified based on macroscopic and microscopic characteristics using stereo microscope (Magnus MSZ-Bi, Noida, India), compound microscope (BX53, Olympus, Tokyo, Japan), confocal microscope (Leica STELLARIS 5, Mannheim, Germany), and scanning electron microscope (S-3400 N, Hitachi, Tokyo, Japan). Furthermore, water-mounted and lactophenol cotton bluestained slides of perithecia, asci, and ascospores were observed and photographed.

. . Scanning electron microscopy (SEM)
Surface morphological studies of O. indica were conducted using SEM. In brief, a small part of the fruiting body was cut with a fine sterile blade, washed with distilled water, and placed on the double-sided carbon tape mounted on the aluminum stub. The mounted sample was left at room temperature, in dust-free condition for drying. Then, this dried sample was coated with gold using a sputter coating unit (E1010 Ion Sputter Hitachi, Tokyo, Japan) at 10 Pa vacuum for 10 s. The images were captured in SEM mode at desired magnifications.

. . Library preparation, sequencing, and phylogenetic analyses
High-quality genomic DNA was isolated from ∼40 mg stroma of O. indica specimen using Wizard Genomic DNA Purification Kit (Promega, Madison, USA) and checked on 0.8% agarose gel. ITS, nrSSU, nrLSU, RPB1, 3 ′ end of TEF-1α, and MAT 1-2-1 were, respectively, amplified with primer pairs described by Chen et al. (2013) and sequenced on a 3730xl DNA Analyzer (Applied Biosystems, CA, USA) via ABI Big-Dye version 3.1 sequencing kit (Applied Biosystems, CA, USA) as per the manufacturer's instructions. The resulting chromatograms were evaluated using Chromas Lite, and a homology search was done using NCBI-BLAST (http://blast.ncbi.nlm.nih.gov). The sequences used for the alignment were downloaded from NCBI and were aligned with Clustal W (Thompson et al., 2003). Pairwise distance matrices and phylogenetic consensus trees were analyzed using MEGA11 software (Molecular Evolution Genetic Analysis; Tamura et al., 2021). The taxon information and GenBank accession numbers used in the molecular analyses are listed in Supplementary Table 1.
Forty-five ITS sequences (Supplementary Table 1A) of Ophiocordyceps species, closely related to O. indica, were used to align and construct a neighbor-joining (NJ) tree (Tamura et al., 2021). NJ estimation of the phylogeny of ITS was performed with 1,000 bootstrap replicates by MEGA11. Beauveria bassiana was selected as outgroup taxa. The Tamura-3-parameter+Gamma distribution with invariant sites nucleotide substitution model was used for the substitution model (Chen et al., 2013). NJ was also constructed for MAT 1-2-1 gene using Kimura 2 parameter nucleotide substitution model. Furthermore, for the construction of a multigene tree, each gene region was independently aligned and improved manually and nrSSU, nrLSU, RPB1, and 3 ′ end of TEF-1α gene sequences were combined to form a dataset. The four gene datasets from the Ophiocordyceps species and datasets obtained from GenBank were then aligned using MUSCLE in MEGA11 (Tamura et al., 2021). Alignments were manually adjusted to allow maximum sequence similarity. Gaps were treated as missing data. Unweighted maximum parsimony (MP) and maximum likelihood (ML) analyses were performed using PAUP * 4.0a10 (Swofford, 2002). The heuristic search option with TBR branch switching and 1,000 random sequence additions was used to infer trees. Maxtrees were 5,000 branches of zero length collapsed, and all multiple parsimonious trees were saved. Clade stability of the trees resulting from the parsimony analyses was assessed by bootstrap analysis with 1,000 replicates, each with 10 replicates of random stepwise addition of taxa (Felsenstein, 1985). Trees were saved and exported to graphics programs. Tolypocladium inflatum was selected as outgroup taxa. Clades that were supported with 50% or greater values were considered significant, while the values below 50% were not displayed in the phylogenetic tree. Furthermore, DNA libraries were constructed with PacBio SMRTbell template preparation kit. SMRTbell template library was checked with Bioanalyzer DNA 12,000 chip (Agilent Technologies, Santa Clara, CA, USA). PacBio Single Molecule Real-Time DNA sequencing technology developed by Pacific Biosciences was considered with 51X coverage for the sequencing of the O. indica genome. The assembly of the long reads obtained was tested on various assemblers and compared to obtain the best assembly of the genome.
. . Isolation of RNA, sequencing, and assembly for transcriptome studies Fresh fruiting bodies of O. indica collected from two different locations were used for the transcriptome study. After collection, these samples were immediately frozen in liquid nitrogen and stored in the laboratory at -80 until further processing. One hundred milligrams of O. indica sample was ground into a fine powder in liquid nitrogen using a mortar and pestle and processed using IRIS method (Ghawana et al., 2011). The purity and concentration of RNA were assessed by determining the absorbance of the sample at 260 and 280 nm using a spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and Agilent Bioanalyzer Chip 7500 Series II (Agilent Technologies, Santa Clara, CA, USA). The RNA was resolved on the gel. The RNA sequencing library was prepared with mRNA-seq 8 sample preparation kit (Illumina, San Diego, CA, USA) using random primers. The concentration ≥800 ng/µl of RNA was used for the preparation of the cDNA library and was run on NovaSeq 6000 Sequencing System (Illumina, San Diego, CA, USA).

. . Assembly and functional annotation
The O. indica fungal genome was first assembled from the most established long-read technology, Pacific Biosciences (PacBio, CA, USA), by using different assemblers: Celera (Denisov et al., 2008;PbcR) and Canu. The best assembler, Canu pipeline along with SSPACE (Boetzer et al., 2011) was used to do an independent whole-genome assembly of raw reads using default parameters to generate the scaffolds. Finally, the completeness of O. indica genome assembly was assessed by the presence of BUSCO orthologous genes [Benchmarking Universal Single-Copy Orthologs; (http://gitlab.com/ezlab/busco)], with the default values. BUSCO uses a set of universal singlecopy orthologs. The assembled transcriptome was then aligned with assembled genome using GMap. The protein-coding genes in O. indica were generated by means of the fully automated Fungal Genome Annotation Pipeline (FunGap; Finn et al., 2016). FunGap predicted gene models in three major steps: (a) pre-processing, (b) gene prediction, and (c) evaluation and filtration. To attain high-quality gene models, the program FunGAP: Fungal Genome Annotation Pipeline (https://github. com/CompSynBioLab-KoreaUniv/FunGAP) was used, where it takes two inputs: O. indica genome assembly in FASTA format and mRNA sequencing reads of O. indica in FASTQ format. This program incorporates multiple gene predictors, viz., Augustus, Braker, and Maker (Stanke et al., 2006;Brandi et al., 2008;Hoff et al., 2016) that assessed all the computationally predicted genes and assemble the best gene model for those particular genes which have shown highest homology to known sequences.
The best score for all predicted gene models of O. indica fungus from FunGAP was then selected on the basis of the alignment of translated protein sequences with well-known functional tools, viz., Pfam, BUSCO, and BLAST (Boratyn et al., 2013). The alignment of the transcriptome, Ascomycota BUSCOs, and homologous peptide to our gene predictions was carried out to evaluate the quality assessment of predicted gene models. InterProScan version 5.16.55 . /fmicb. . (Zdobnov and Apweiler, 2001) was then used to analyze the functional annotation of the predicted gene models. Non-coding RNA (ncRNA), including transfer RNA (tRNA), ribosomal RNA (rRNA), small nuclear RNA (snRNA), and small nucleolar RNAs (snoRNA), were predicted using cmsearch in Infernal using the Rfam library Rfam.cm.1_1 with the hits selected above the evalue threshold of 1e-5.

. . Metabolite profiling
Different extraction methods viz., cold water, hot water, and methanol (50:50, 80:20, and 100%) were tried for the extraction of O. indica and O. sinensis samples using an ultrasonic water bath (Elma S 300 H, Elma Schmidbauer GmbH, Singen, Germany). Metabolomic studies were done using UHPLC-Q-ToF-IMS (Agilent, Santa Clara, USA) with an Eclipse Plus C18 RRHD column . The identification of unknown metabolites other than targeted ones was searched using acquired data against the METLIN database. Identified metabolites were classified according to their biochemical groups namely, sugars, fatty acids, amino acids and their derivatives, nucleosides, and others. Amino acids (Agrawal et al., 2019) and mannitol were quantified in hot water extract using UPLC.

. . Antioxidant activities
DPPH and ABTS radical-scavenging activities of hot water extracts of O. sinensis and O. indica were performed (Blois, 1958;Joshi et al., 2015Joshi et al., , 2018. The reaction was performed using different concentrations of each sample (1-10 µl) by adding 2.9 mL of 0.06 mM DPPH solutions to the reaction mixture and vortexed. All the tubes were covered with aluminum foil to carry out the reaction in the dark. The reaction mixture was incubated for half an hour at 25 • C, and the absorbance was recorded at 517 nm wavelength. Similarly, the ABTS assay was performed using an aqueous solution of 2,2 ′ -azinobis-(3-ethylbenzothiazoline-6-sulphonate; ABTS, 7 mM, 5 mL) and potassium persulphate solution (88 µL, 140 mM) which were mixed and kept overnight in the dark to generate ABTS .+ cation (2.45 mM). ABTS .+ solution is diluted with ethanol to achieve an absorbance of . /fmicb. . 0.7 ± 0.05 at 734 nm. The reaction was performed using different concentrations of each sample (1-10 µl) which were mixed with 2 mL of diluted ABTS .+ radical solution, and the absorbance was recorded at 734 nm after 6 m of initial mixing. The experiments were carried out in triplicates, and IC 50 values were calculated by plotting a graph between absorbance and concentration.

. . Cell viability studies
The cellular viability of peritoneal macrophages in response to hot and cold water extracts of O. indica was assessed using four different concentrations (25, 50, 100, and 200 µg/mL;Sharma et al., 2018).
The viability of cells was determined using the following equation: where test sample = OD of cells treated with tested samples and control sample = OD of cells treated with only media.

. . . Nitrite production
The peritoneal macrophages (0.5 × 10 6 ) were supplemented with four non-lethal concentrations (25, 50, 100, and 200 µg/mL) of hot and cold water extracts O. indica and stimulated with LPS (1.5 µg/mL; Sigma Aldrich, St. Louis, USA) for 16 h in a CO 2 incubator. The supernatant was collected after incubation, and NO production was analyzed using Griess reagent as per the manufacturer's protocol (Promega, USA; Sharma et al., 2018).

. . Statistical analyses
The data were statistically analyzed using the SPSS software (SPSS version 27.0.1). Student's ttest was used to determine and compare the means (McDonald, 2014 Sexual morph: Stroma solitary, simple, or branched (up to 120 mm long above soil, dia 0.42-1.70 mm and 6-24 mm below soil), protruding from head between two eyes of host caterpillar (20-36 mm; Figures  . . Genome sequencing, assembly, and evaluation ITS has been used as an important molecular marker in the identification and genetic analysis of fungi including Ophiocordyceps species. In the present study, molecular data and systematic analyses based on the ITS region of the Ophiocordyceps specimen (KX679571) revealed its 97.67% similarity with Hirsutella uncinata ( Figure 3C). Furthermore, a multigene phylogeny study based on nrSSU (MZ571406), nrLSU (KY486751), RPB1 (MZ514912), 3 ′ end TEF-1α (MZ514914), and MAT 1-2-1 (MZ514913) revealed it to be a novel species of the genus Ophiocordyceps ( Figures 3D, E). The four-gene combined dataset consisting of 2,828 base pairs (nrSSU 812 bp, nrLSU 635 bp, TEF-1α 911 bp, and RPB1 470 bp) of O. indica were used to build a multigene phylogenetic tree of closely related species of Ophiocordyceps ( Figure 3D; Supplementary  ). The data generated on the PacBio were long reads with an average length of 5,807 bp with ∼6.9 Gb in overall 1,182,152 in fasta read sequences with minimum and maximum read lengths of 35 and 51,733 bp, respectively (Accession No. PRJNA699552). Canu assembler (Koren et al., 2017) was found to be the best in terms of the lowest number and size of contigs (356 and 59,047,143 bp) with an N50 contig length of 6,27,704 bp. Finally, 229 scaffolds (13,118 SSRs) were generated after scaffolding with a total size estimated to be 59,247,067 bp and a ScafN50 value of 852,095 bp (Table 1). BUSCO (Simão et al., 2015) [benchmarking universal single-copy ortholog] evaluation of the O. indica genome revealed that 95.8 and 3.3% of the 1,315 expected Ascomycota BUSCO conserved genes were identified as complete and fragmented, respectively. A few numbers of missing genes (0.9%) indicated the completeness of the O. indica genome ( Table 2). Out of these 13,118 SSRs, the most abundant type was mono-nucleotides repeats (5,255; 40.05%), followed by di-nucleotides (4,022; 30.7%), tri-nucleotide (3,328; 25.36%), tetra-nucleotides (191; 1.45%), hexa-nucleotides (188; 1.43%), and penta-nucleotides (134; 1.02%' Table 3 and Figure 4A).
To annotate gene models, an average of ∼22 Gb of clean RNAseq (sequencing) data were generated from the Illumina Nova-Seq sequencing platform and obtained 25,558 transcripts (average length 2,096 bp). The assembled transcriptome showed more than 99.46% alignment with a criterion of using >90% in terms of coverage and identity, again indicating a very good quality of the assembled genome. The annotation software, FunGAP (Min et al., 2017), used multiple gene predictor models (Stanke et al., 2006;Brandi et al., 2008;Hoff et al., 2016) along with transcripts as a query finally resulted in 10,613 protein-coding genes. Using BUSCO from Ascomycota lineage, the best-annotated protein models were observed with (1,264 out of 1,315) 96.10% as complete, (51) 3.90% as duplicated, (49) 3.72% as fragmented, and (2) 0.15% as missing, revealing a very good quality of our gene predicted models ( Table 2). The gene prediction in O. indica produced by FunGap revealed a higher quality and hence can be processed for downstream protein analysis. InterProScan (Quevillon et al., 2005) findings revealed that 69.2% of annotated models had matches in the InterPro database, 6.7% in KEGG, and 49.1% in Gene Ontology (GO ; Table 4). Non-coding RNA (ncRNA) genes using Infernal (Nawrocki and Eddy, 2013) Figure 4C). Pfam annotations identified 10,718 Pfam domains within 10,613 genes in O. indica, which were within the same range (∼1% of annotated Pfam) as O. sinensis (9,343 domains within 7,939 genes). Furthermore, the numbers of unrepeated Pfam entries were almost similar such as 3,874 (36.5%) for O. indica and 3,339 (42.05%) for O. sinensis. The most frequent Pfam domains common in both genomes were "WD domain, G-beta repeat, " "Protein kinase domain, " "Major facilitator superfamily, " "Mitochondrial carrier protein, " and "Fungal Zn(2)-Cys(6) binuclear cluster domain." These observations suggested that they share most of the functions that can determine their lifestyles. The abovementioned Pfam domains represented more than 1% of O. indica fungal genome; however, only two Pfam domains ("WD domain, G-beta repeat" and "Protein kinase domain") were reported to have more than 1% in O. sinensis.
In contrast to lower number (7,939) of protein-coding genes (Xia et al., 2017)   similarity among Ophiocordyceps species ( Figure 4E). Nearly 2,684 and 496 genes specifically belong to O. indica and O. sinensis, respectively, wherein 1,918 and 102 proteins were annotated with functional terms based on InterProScan analysis. Among the specific genes in O. indica, the most dominant Gene Ontology terms were associated with oxidation-reduction process (GO: 0055114), transmembrane transport process (GO: 0055085), pathogenesis (GO: 0090729), and DNA-dependent regulation of transcription (GO: 0043565). The finding of oxidation-reduction processes and the trademark for parasite-host interactions including cytochrome P450, FAD, and NAD-binding domains were specific and over-represented in the O. indica genome. Other protein-coding genes, viz., membrane transporters/branched amino acid metabolism (PF00528/PF02653), Zn (2)-C6 fungaltype (PF00172), LysR bacterial regulatory protein (PF00126), and heat-labile enterotoxin alpha chain domain (PF01375), were most exclusively present.
On the contrary, five genes of a DNA-binding motif (NUMOD3 motif; 2 copies; PF07460) present in the genome of O. sinensis were found to be totally absent in the O. indica genome. Similarly, a total of 328 syntenic blocks were identified between O. indica and C. militaris comprising 9,265 collinear genes (∼46%) which were also in the same range, considering collinearity genes of O. sinensis and C. militaris.

. . Comparison of pathogenesis proteins in O. indica fungus with other entomopathogenic fungi
Based on peroxidase profiles, we explored the peroxidase genes amid species of interest using "hmmsearch" program executed in Hmmer (version 3.1b1) package with default parameters. As a result, a total of 54 (0.    Figure 4F), proposing that a higher number of peroxidase genes assist in ROS detoxification in both Ophiocordyceps species. It is well-reported that secreted enzymes are upregulated during host-pathogen interaction (de Bekker et al., 2017), which includes many pathogenicityrelated genes (proteases and heat-labile enterotoxin alpha chain domain). Using MEROPS, a total of 343 (3.23%) protease genes were identified in the O. indica (Supplementary Table 2A). In comparison with other fungal species, a comparable number of membrane transport proteins (399) were found in O. indica (Supplementary Table 2B). Using the KinBase database, a total of 75 kinases were identified in O. indica (Supplementary Table 2C).
The heat-labile enterotoxin alpha chain domain (PF01375) was highly amplified in all entomopathogenic fungi of the family Ophiocordycipitaceae except O. sinensis ( Figure 4F). A large repertoire of 29 (0.27%) heat-labile enterotoxin alpha chain domains was found in the O. indica genome in comparison with 9 (0.11%) in the O. sinensis ( Figure 4F). Enterotoxins are considered important in host-pathogen interactions (Dong and Yao, 2012). Fifteen proteins from a total of 29 encoding heat-labile enterotoxins (alpha chain) were found specific only in the O. indica genome and not exclusively shared with O. sinensis, reflecting the specificity of the presence of a higher number of enterotoxins, playing an important role in manipulating the behavior of host insects (Tudzynski et al., 2012).
Many insect pathogens require carbohydrate-active enzymes (CAZymes), which are liable to infect their hosts by degrading the cuticle (Cantarel et al., 2009

. . Orthology assignment and phylogenetic analysis
The OrthoFinder analysis of 1,95,218 genes from 19 fungal species (O. indica, O. sinensis, O. australis, O. unilateralis, O. polyrhachis-furcata, O. camponoti-rufipedis, H. sinensis, T. inflatum, M. anisopliae, M. acridum, C. militaris, and B. bassiana as 12 insect-fungal species; F. graminearum, V. alfalfa, M. grisea, G. clavigera, S. sclerotiorum, and B. cinerea as 6 plant-fungal species; and Saccharomyces cerevisiae as outgroup) identified 13,076 gene families/orthologous groups. The total number of genes in the orthologous group was found to be 173,334 (88.8%) while failing to allocate 21,884 (unique) protein-coding genes into any group. One thousand two hundred and thirty-three orthologous groups consisted of a single-copy protein from all fungal species, while an additional 2,234 orthologous groups containing multiple members in one or more species were considered for phylogeny, using RAxML along with S. cerevisiae as an outgroup species ( Figure 4G). The phylogenetic analysis revealed that the O. indica can be considered as the novel species in the family of Ophiocordycipitaceae with existing members O. sinensis, O. unilateralis, O. australis, O. polyrhachis-furcata, O. camponoti-rufipedis, H. sinensis, and T. inflatum (Figure 4H).

. . Evolution of expanded and contracted gene families in O. indica
Computational Analysis of gene Family Evolution (CAFE; Wang and Wang, 2017) identified 117 gene families (with 361 genes) in O. indica with significantly higher than expected rate of gains/losses (P ≤ 0.05) among 14 fungal species ( Figure 4I). Ophiocordyceps indica exhibited expansion (significant P < 0.01) in gene families, responsible for fungal pathogenicity that comprised of peroxidase (PF01328), heat-labile enterotoxin alpha chain (PF01375), lipase class 3 (PF01764), short-chain dehydrogenase (SDR; PF00106), and N-terminal domain on NACHT_NTPase and P-loop NTPases (PF17107). The highly enriched peroxidase genes were also found to be expanded in O. sinensis (Xia et al., 2017); this reflected a strong need for peroxidase genes in both Ophiocordyceps species for adapting to the prevailing Himalayan environments. However, those over-represented pathogenic proteins in O. indica, such as heat-labile enterotoxin alpha chain and lipase class 3 proteins, were not found to be expanded in O. sinensis.

. . Metabolite profiling
Ophiocordyceps species possess different bioactive constituents, viz. nucleosides, amino acids, polysaccharides, sterols, proteins, and polypeptides, which make this genus a fungus of medicinal and nutraceutical importance. Therefore, we targeted eight major marker constituents (adenine, adenosine, cordycepin, guanosine, inosine, thymidine, thymine, and uracil) using UPLC-QTOF-MS (Table 5; Supplementary Figures 1-4). The compounds were identified by comparing their retention times and UV spectra with those obtained on injecting standards under the same conditions or by spiking extracts with stock standard solutions. Amazingly, all eight nucleosides were present in O. indica and O. sinensis. The hot water extraction (2,303.25 and 2,579.83 µg/g) was more efficient than methanolic (80:20) extraction (393.42 and 482.23 µg/g) in O. indica and O. sinensis, respectively (Table 5) (Table 5). In both species, guanosine was found to be the maximum and thymine in the lowest quantities. Thymine was below the quantification level in the hot water extract of O. sinensis and was not detected in alcoholic extracts in both O. indica and O. sinensis.
A comparable quantity of amino acids (6.15%) was quantified in O. indica and is comparable to O. sinensis (7.77%). Out of the targeted 14 amino acids, the most abundant were norvaline, arginine, leucine, proline, and valine in O. indica (Table 6)  Values are expressed as mean ± standard deviation (SD; n = 3). The same letters in the same column of each, hot water and methanolic extract are not statistically significant according to Student's t-test with a P-value of ≤ 0.001. Ala, alanine; Arg, arginine; Asp, asparagine; His, histidine; Leu, leucine; Lys, lysine; Met, methionine; Nva, norvaline; Pro, proline; Ser, serine; Thr, threonine; Trp, tryptophan; Tyr, tyrosine; Val, valine. Values are expressed as mean ± standard deviation (SD; n = 3). The same letters in the same column are not statistically significant according to Student's t-test with a P-value of ≤0.001.  Values are expressed as mean ± standard deviation (SD; n = 3). The same letters in the same column are not statistically significant according to Student's t-test with a P-value of ≤0.05.
nucleosides as the major classes of non-targeted metabolites in the extracts of O. indica. The abundance of identified metabolites was visualized using heat maps, illustrating the metabolite variations (Supplementary Figure 5).

. . Biological activity and safety studies
Ophiocordyceps indica exhibited promising antioxidant activity suggesting its therapeutic potential. These activities were statistically at par with O. sinensis (Table 7). The IC 50 DPPH was 1.21 times higher in O. indica (2.77 mg/mL) as compared to O. sinensis (3.34 mg/mL); however, the reverse trend was observed in IC 50 ABTS, where O. sinensis (0.88 mg/mL) was 1.05 times more efficacious than O. indica (0.92 mg/mL). The present findings are in accordance with earlier findings (Singh et al., 2018).
In vitro study revealed that hot water extract of O. indica was safe and non-toxic and did not compromise cell viability ( Figure 5). The extracts of O. indica enhanced the nitric oxide (NO) production which is considered responsible for many biological activities including aphrodisiac activities, by virtue of which O. sinensis is known as Himalayan Viagra. LPS-stimulated macrophages significantly (p ≤ 0.05) enhanced NO production (13.24 ± 0.53 µM) in comparison with control (vehicle control; 4.34 ± 0.34 µM), while hot water extract resulted in higher NO production as compared to cold water extract. These results established the potential of this novel O. indica as an activator of NO production.

. Discussion
We reported a novel species, O. indica, from the treeline areas (2,202-2,653 m AMSL) of H.P., Indian Western Himalayas. Many Ophiocordyceps species have only slight morphological differences; therefore, molecular phylogenetics plays a significant role in confirming their status as separate evolutionary species (Kobmoo et al., 2012;Araújo et al., 2018;Tasanathai et al., 2019). The phylogenetic trees using NJ, ML, and MP methods suggested that O. indica is a close relative of O. xuefengensis, O. sinensis, O. lanpingensis, O. stylophora, and H. uncinata. The multigene analyses also confirmed its distinctiveness among other Ophiocordyceps spp. Generated morphological, genetic, chemical, and biological activity data showed its resemblance with the alpine meadow (3,000-5,500 m AMSL) medicinal fungus, O. sinensis. Based on habit, habitat, morphology, ITS, multigene (nrSSU, nrLSU, RPB1, 3 ′ end TEF-1α, and MAT 1-2-1) phylogenetic, and whole-genome sequence studies, it is identified as a novel species of the genus Ophiocordyceps, for which the name Ophiocordyceps indica sp. nov. is proposed. Morphologically, this fungus appears somewhat like O. lanpingensis, O. xuefengensis, O. robertsii, and O. stylophora; however, it is distinct from these species (Table 8). The stroma of O. indica is much smaller than O. xuefengensis, whereas O. stylophora is reported to infect coleopteran larvae and we have recorded it from a lepidopteran larva, Thitarodes sp. (Accession No. OR044664). In O. robertsii, the fertile part is present at the top of the stroma, whereas in O. indica it is present in the middle. Ophiocordyceps indica looks quite similar to O. lanpingensis, however, it differs in having a hardy stroma as compared to wiry to pliant or fibrous in later (Chen et al., 2013). In the present study, the perithecium of O. indica is smaller than O. robertsii and O. sinensis (Chen et al., 2013). Asci were filiform in O. indica as compared to cylindrical in O. xuefengensis (Wen et al., 2013). Ascospores were thread-like, multiseptate, hyaline, and filiform (68.0-307.7 × 3.03-5.54 µm) in O. indica and are larger in size as compared to filiform ascospores (34.5-48.0 × 2.0-2.4 µm) of O. acicularis (Liang, 2007). Filiform, multiseptate ascospores of varying sizes (8.0-470.0 × 1.5-10.0 µm) have been reported in O. sinensis (Shrestha et al., 2010). Comparative characteristics of O. indica and seven phylogenetically closely related species are presented in Table 8, which clearly depicts its distinctiveness from other spp. Thus, the identification and introduction of O. indica (a sister taxon of O. sinensis), a novel species from low-height areas, are of great significance as it can be an alternative to highly prized O. sinensis.
Ophiocordyceps sinensis is one of the world's most expensive natural medicinal resources, which has and is providing bread and butter for millions of traditional collectors in Himalayan regions (Li Y. et al., 2021;Smith-Hall and Bennike, 2022). Nevertheless, due to habitat degradation, climate change, overexploitation, grazing, and anthropogenic pressures, this fungus has reached an alarming scale of extinction, resulting in it being listed as vulnerable on the IUCN Red List of Threatened Species and the Red List of China's Biodiversity-Macrofungi (Yao et al., 2020;Wei et al., 2021). One of the most common reasons for the reduction in the yield of insect-fungal complex is that the harvesters almost collect all the Ophiocordyceps that they encounter, which might result in the production of fewer spores for infecting new caterpillars in the subsequent years. To the best of our knowledge, O. sinensis has never been reported from any of the geographical locations, including the alpine meadows and treeline areas of H.P., India, despite resemblances of topography and climatic conditions to that of other locations of its distribution within and outside India. Therefore, the identification of this novel O. indica is of tremendous importance and can be exploited for the upliftment of the socioeconomic status and livelihood of the rural inhabitants of H.P. through scientific interventions.
If we compare the habitat, O. indica is recorded from tree line areas from lower heights, in contrast to the alpine meadow habitat of O. sinensis, while the distribution of O. lanpingensis is reported from mountain grasses (2,000-3,400 m AMSL) over Hengduan mountain region, China. In the present study, trees (Picea smithiana, Acer palmatum, Aesculus indica, and Prunus cornuta) and lower plants (Carpesium sp., Diplazium sp., Dryopteris . /fmicb. . indica. Results are expressed as mean ± SD (standard deviation), where n = . The statistical significance of the data is presented as compared to the control; ns represents that a non-significant di erence was observed. * and **** represent statistical significance with P < . , and P < . , respectively.
The immature of O. indica at the time of collection, as well as after drying, was morphologically quite similar to O. sinensis except the color of the caterpillar part and the thickness of the stroma. Caterpillar was little dark brown in O. indica as compared to the golden color of O. sinensis and the lighter shade of O. lanpingensis (Chen et al., 2013). The variation in color from light brown to brown or brownish-black may be attributed to habitat differences as well as melanin concentration due to environmental factors (Dong and Yao, 2012;Shrestha et al., 2012). During the study, we were able to collect a live lepidopteran caterpillar from the soil that was light brownish in color similar to the ghost moth caterpillar, the host of O. sinensis as reported in the literature ( Figure 1F).
The size, color, weight, stiffness, smell, taste, and robustness of intact caterpillar and stroma decide and influence the demand, market price, and public acceptability of the Ophiocordyceps (Zhou et al., 2014;Li Y. et al., 2021), where immature individuals are more valuable (Hopping et al., 2018). Golden-colored and vigoroussize pieces are in great demand in the international market, and generally, 2,000 pieces/kg is considered the best grade (A++) Ophiocordyceps (Supplementary Table 4). The best, healthy, and intact specimens of O. indica were sorted, weighed, and counted, and excitingly, O. indica samples (1,661 pieces/kg) qualify for the best category for fetching high market price.
Our next curiosity and challenge was to find the availability and abundance of O. indica; therefore, surveys were conducted in subsequent years and information was gathered from some local inhabitants by showing the collected specimens. Based on the data and information gathered from a few locations only, it is inferred  that there is a huge quantity (in Kgs) of O. indica in the nearby unexplored areas that was not even collected due to ignorance regarding this highly prized fungus and unavailability of the local market, in addition to fear of administrative rules and regulations for its collection. As of now, there are no guidelines for collection of this novel O. indica as it is not reported and recorded earlier. Therefore, it is necessary to frame rules, regulations, and policies, especially for sustainable harvesting and conservation in its natural habitat to avoid overexploitation and overharvesting, as we all witnessed such issues in the case of O. sinensis. This can be achieved by declaring some protected areas where harvesting can be totally prohibited or restricted.
Ophiocordyceps species are known to be parasitic on many species of insects with intriguing host selectivity. Recently, it was suggested that the parasitic association could have resulted in greater selective pressure to adapt the host's immune response in highly specific niches. The phylogenetic trees suggest that Ophiocordyceps species are not a monophyletic group. Furthermore, the analyses of the O. indica genome and transcriptome unveiled that gene families associated with fungal pathogenicity have experienced rapid evolution and undergone natural selective pressures. In the present study, of the 22,258 transcripts assembled, ∼25,420 (∼99.46%) had matches with our genome assembly having parameters >90% coverage and identity. A higher percentage provided by the benchmark approach, BUSCO, in terms of both transcriptome and genomes, validated a high quality of the genome assembly, securing further comparison with entomopathogenic fungi including O. sinensis. In contrast to O. indica, a global contraction of gene families was apparently observed which evidenced the removal of non-collinear genes and a total loss of gene families in O. sinensis.
Insect hosts often expeditiously produce an assload of reactive oxygen species (ROS) against the pathogens, and during evolution, a ROS antioxidant defense system is developed in which peroxidases (ROS-scavenging enzymes) are well-known and fundamental components (Tudzynski et al., 2012). In addition, carbohydrate-active enzymes (CAZymes), responsible for the synthesis, degradation, and modification of all carbohydrates .
/fmicb. . present in the genome including glycoproteins and glycolipids, have also been identified in our study. The functional properties of these gene families revealing expansion status are enriched in the transport process (ABC transporters) and energy metabolism that might have been lost in O. sinensis due to accumulation at higher altitudes for the need of conservation of energy. The O. indica genome also revealed substantial expansion of gene families that are known to participate in fungal pathogenicity, including peroxidase activity, cytochrome P450, and, especially, heat-labile enterotoxin that was not found to be expanded in the O. sinensis genome. The heat-labile enterotoxin (alpha chain) is known to exhibit manipulation of host insects by interfering with the production of chemical signaling molecules and inferred specificity of the O. indica genome. Lipase proteins are involved in the degradation/detoxification of cuticular lipids of insect hosts, thus, inferring their higher importance for pathogenicity (Wang and Wang, 2017). Additionally, these can reduce the energy required to carry out catalytic activity reactions at unusual temperatures. Fungal melanin biosynthesis is dominated by the enzymes involving short-chain dehydrogenase (SDR), which are required for pathogenicity in M. oryzae (Wang and Wang, 2017), and thus may play a pivotal role in O. indica fungal pathogenesis.
To counter harsh environmental conditions, especially UV radiations, fungi will require UV protectants and heat stabilizers, such as the tyrosinase gene for normal functioning of the cell. A considerable expansion in tyrosinase-enriched gene families (PF00264) in the O. indica genome might help in providing an effective way for the stability of spores and resistance to UV radiations (Shang et al., 2012). Membrane transporters are implicated not only in the overall transport of the pathogenic and virulence compounds but also in protection against insect defense compounds as secondary metabolites, probably by exporting the host-derived antimicrobial compounds outside the fungal cell (Merzendorfer, 2014). As expected, highly expanded genes were found to be involved in the transportation process that included ABC transporters (PF00005) and binding-protein-dependent transport system inner membrane component (PF00528) in the genome of O. indica. This finding was in contrast to O. sinensis where the functional properties of contracted gene families were enriched with transport processes and energy metabolism (Xia et al., 2017).
Some of the enriched PFAM functions were also found related to DNA transposition (Retrotransposon gag protein; PF03732), protein-protein interaction (ankyrin repeats; PF12796), and carbohydrate-binding domain (GLEYA adhesin domain; PF10528). In contrast, gene families exhibiting contraction status (P < 0.01) were mainly involved in the pathogenesis such as the WD domain, G-beta repeat (PF00400), and NACHT (PF05729). This finding was in contrast to the expansion of both Pfam domains in O. sinensis (Xia et al., 2017); however, the loss of both genes can be thought of compensated by the higher expansion of other pathogenic proteins (heat-labile enterotoxin alpha chain) in the genome of O. indica.
In addition, in fungi, mating and sexual development regulators are encoded by mating-type genes. Different sets of mating-type genes are requisite by heterothallic ascomycete species to control non-self-recognition and mating of compatible partners of different mating types. Homothallic (self-fertile) species also carry matingtype genes in their genome that are essential for sexual development (Klix et al., 2010). Generally, the mating system is controlled by the mating-type (MAT) locus in Ascomycota fungi (Klix et al., 2010). MAT1-1-1 and MAT1-1-3 play a role in fruiting body development, whereas MAT1-1-2 is essential for sexual reproduction (Klix et al., 2010). O. indica possessed MAT1-1 idiomorph including MAT1-1-3/A-3 mating-type gene and MAT1-2 idiomorph having MAT1-2-1 mating-type gene, concluding it to be a homothallic species, the second species of its own type after O. sinensis in the family Ophiocordycipitaceae.
In addition to identifying any new Ophiocordyceps spp., the identification of marker compounds is the utmost requirement for establishing their economic importance. Therefore, metabolomic studies of O. indica were conducted and compared with O. sinensis. Hot water extraction was the best extraction method and yielded 5.85 times more nucleosides as compared to alcoholic extraction in O. indica (Table 5), while it was 5.35 times higher in O. sinensis. These results suggest that O. indica can be consumed in tea, soup, hot drinks, or other cooked items, just like O. sinensis (Chen et al., 2010;Singh et al., 2018;Kaushik et al., 2019). Mannitol, adenosine, and cordycepin are the important qualitative parameters, and our results showed that all these were present in both species. Most of the targeted metabolites were present in both O. indica and O. sinensis; however, their content varied, possibly due to the difference in metabolism and transformation for nucleosides and amino acids (Chen et al., 2018). Cordycepin is also reported from other species, namely O. lanpingensis and O. xuefengensis (Chen et al., 2013;Jin et al., 2018). Mannitol was (10.03%) 1.12 times higher in O. indica than O. sinensis and is almost equal to that reported (10.04%) in natural O. sinensis of Sichuan Province, China. Cordycepin content (0.06%) in O. indica is the same, while adenosine (0.24%) is 4 times higher than previous records on O. sinensis (Zhou et al., 2019). In O. indica, adenosine (243.32 µg/g) is higher than in O. lanpingensis (113.11 µg/g; Chen et al., 2018). In both O. indica and O. sinensis, guanosine was maximum as reported earlier in natural and artificial C. sinensis (Li et al., 2001). The presence of guanosine has also been reported from O. lanpingensis (Chen et al., 2013). Our targeted metabolite approach is further supported by METLIN results, which revealed the presence and matching of 80% of compounds of both species (Supplementary Figure 5). If we look at O. indica from a chemical composition point of view (targeted and non-targeted, i.e., METLIN search approach), the presence of a comparable amount of nucleosides, amino acids, mannitol, and other chemical constituents as compared to O. sinensis suggests that it can be utilized as an acceptable dietary supplement or pharmacological product after generating necessary safety data.
To the best of our knowledge, this study is the first report of a new species of Ophiocordyceps (O. indica) from the lowheight treeline areas of the Indian Western Himalayan region supported by morphological, whole-genome, and phylogenetic analyses. The results exemplified that O. indica is clustered in the genus Ophiocordyceps of Ophiocordycipitaceae and is closely related to O. sinensis and H. sinensis. The study presents new insights into the taxonomy, genetics, systematics, and evolution of this important O. indica. The comparative analysis with genomes of Ophiocordyceps spp. and other fungi will enrich our understanding The identification of this novel O. indica from low-height approachable sites will open new avenues for studying biology, host insect-fungus interaction (host screening, breeding, infection rate, and fruiting body development mechanism), cultivation, and other ecological aspects that are difficult to study in other Ophiocordyceps spp., especially O. sinensis due to its habitat in tough terrains. In addition, many more aspects of this species need to be researched, which will definitely provide new research outcomes. Such indepth studies will undoubtedly strengthen our findings and lead to further expansion of the subject to exploit O. indica as a potential alternative to highly prized medicinal mushrooms, O. sinensis "Himalayan Viagra". These endeavors will help the rural population of the region, by generating employment, thereby improving their livelihood in addition to conserving foreign cash reserves and generating revenue.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

Ethics statement
In the present study, sample collection was conducted according to India's Biological Diversity Act 2002. This act permits its bonafide citizen to use biological resources for scientific research (Venkataraman, 2009).

Author contributions
AS: wet lab experiments, data analyses, writing original draft, and manuscript editing. EK: bioinformatics analysis. RJ: chemical profiling and metabolomic studies. PK: experiments. AK: web link and data arrangement on website. DK: chemical profiling, metabolomics, and data analysis. MS: genome sequencing. VA: conceptualization, planning, computational analyses, and writing. GN: conceptualization, planning, execution of the survey, collection of samples, morphological studies, and manuscript editing. All the authors have read, critically reviewed, and approved the submitted version. . /fmicb. .