Diverse Upregulated Transcripts During Conidiation of the Aphid-Obligate Fungal Pathogen Conidiobolus Obscurus (Entomophthoromycotina)


 Filamentous fungi in the order Entomophthorales are the main natural regulators of insect populations. Conidiation is crucial for entomopathogenic fungi to explore host resources due to the multifunction of conidia such as growth, infection, and stress resistance; however, the molecular mechanisms underlying the conidial functions in Entomophthorales is unknown. This study analyzed the differentially expressed transcriptomic patterns in three conidiation stages (pre-conidiation, emerging conidiation, and post-conidiation, respectively) of the aphid-obligate pathogen Conidiobolus obscurus (Entomophthoromycotina). The emerging conidiation stage vs. pre- or post- conidiation stage had 3,091 and 3,235 differentially expressed genes (DEGs), respectively, wherein 2,915 upregulated DEGs were putatively related to the conidial functions. A weighted gene co-expression network analysis showed that 772 hub genes in conidiation, which were related to cuticular component degradation, oxidative phosphorylation, ribosomal biogenesis, cell wall/membrane biosynthesis, MAPK signaling pathway, secondary metabolite biosynthesis, and other metabolic processes. This implied that the conidia of Entomophthorales have abundant transcripts with various functions to favor a quick response to the surrounding environment and effectively explore the host resources.


Introduction
Filamentous fungal species undergo asextual reproduction in their life cycle by the formation and release of asextual spores from aerial hyphae (Jung et al. 2014;Ruger-Herreros and Corrochano 2020). Asexual spores are formed by an individual fungus through mitosis and subsequent cell division (Park and Yu 2012). Conidiospores (conidia) are the main asexual fungal spores that favor fungal spread, survival, and colonization in new habitats (Ball et al. 2019). Pathogenic fungi also use conidia as infectious propagules to recognize and infect host plants and animals (Pell et al. 2001). Thus, fungal conidia have multiple functions and are responsible for adaptation in response to different abiotic and biotic environments.
Entomopathogenic fungi are widely distributed in the environment and regulate host insect populations by causing mycosis and epizootics (Hajek 1997;Pell et al. 2001;Gryganskyi et al. 2017). The entomopathogens in the orders Hypocreales and Entomophthorales develop millions of differently shaped clonal ballistic conidia in their infection cycles to explore host resources, with the capability of adapting to different environments (Spatafora et al. 2016). For example, the conidia of Hypocreales such as Metarhizium anisopliae and Beauveria bassiana undergo a dimorphic transition of hyphae in the host insect integument and blastospores in the hemocoel (Jung et al. 2014). The conidial germination of entomophthoralean fungi depends on the interface they touch. After landing on non-host surfaces, successive conidia are quickly ejected to increase the chances of encountering hosts (Zhou et al. 2014a). Once the host is recognized, the conidia germinate into infectious germ tubes and penetrate the host cuticles for several hours under warm and moist conditions (Grell et al. 2011;Zhou et al. 2014a). After exhausting host nutrition, the fungi re-sporulate to start a new infection cycle. Conidia-based formulations of Hypocreales for insect control have been applied, but the formulation and application of Entomophthorales remain stagnant (Jackson et al. 2010). Furthermore, little attention has been paid to the molecular mechanism of Entomophthorales conidiation.
Variation in gene expression is the basis for phenotypic heterogeneity among mycelia and conidia (Wang et al. 2021). However, the air-borne conidia are generally dormant when they detach from conidiophores. The mechanism of the quick response of dormant conidia landing on host or non-host surfaces of entomopathogens remains unclear. The previous studies showed that host-mimicking media plus host components can induce diverse expression of genes encoding enzymes such as glycoside hydrolases, lipases, and proteases for penetrating the host cuticle (Xu et al. 2006;Lopez-Perez et al. 2012). Thus, pro ling gene expression patterns during the conidiation period favors understanding the infection mechanism and the speci c roles of conidia. This study aimed to screen the conidiation-related genes that are involved in the pre-conidiation (pre-C), emerging conidiation (C), and post-conidiation (post-C) stages of the aphid-obligate fungal pathogen Conidiobolus obscurus (Entomophthoromycotina) using transcriptomic analysis. A weighted gene co-expression network analysis (WGCNA) was used to identify the hub genes for understanding the molecular characteristics of conidiation.

Materials And Methods
Fungal culture and scanning electron microscopy C. obscurus 7217 was obtained from the United States Department of Agriculture-Agricultural Research Service Collection of Entomopathogenic Fungal Cultures (USDA-ARSEF, Ithaca, NY, USA) and stored at −80°C (Zhou et al. 2014b). The isolate was cultured on rich Sabouraud dextrose agar plus yeast extract (40 g L − 1 dextrose, 10 g L − 1 peptone, 10 g L − 1 yeast extract, and 15 g L − 1 agar) for 4 days in petri dishes at 24 ± 1°C with a 12:12 h light: dark photoperiod. Culture pieces were mashed and transferred into 50 mL liquid media (150-mL ask) and incubated in a shaker at 120 rpm and 24 ± 1 °C for 3 d. Fresh mycelia (pre-C stage) were collected using a 0.2-µm lter, and evenly poured into a 90-mm Petri dish to form a mycelial mat while removing any excess water using sterile paper. Moist conditions (100% relative humidity) were maintained for 12 h at 24 °C to form the C stage. The mycelial mats maintained for 72 h at the same temperature and humidity represented the post-C stage.
For observing the morphological characteristics of C. obscurus in different stages of conidiation, the mycelial samples were collected. After xation, dehydration, and platinum ion plating, the specimens of the mycelial mats of three conidiation stages were separately observed with Philips XL30-ESEM environmental scanning electron microscope (SEM) (Zhou et al. 2014a).

RNA extraction and transcript assembly
Three biological replicates were sampled from each of the three stages: pre-C, C, and post-C to screen for conidiation-related genes in C. obscurus. Total RNA was extracted using the RNAiso Plus kit (TaKaRa, Tokyo, Japan), and its concentration was measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scienti c, New York, NY, USA). RNA degradation and contamination (especially DNA contamination) were monitored on 1% agarose gels and samples were sent to Biomarker Technologies Co., Ltd. (Beijing, China) for transcript sequencing. Ribosomal RNA was removed from total RNA using a Ribo-Zero rRNA removal kit (Epicenter, Madison, WI, USA). A total of 1.5 µg rRNA-free RNA per sample was used to generate a sequencing library using the NEBNext Ultra Directional RNA library prep kit for Illumina (New England BioLabs, USA) on an Illumina HiSeq 4000 platform (BGI, Beijing, China) according to the sequencing pipeline ( Figure S1A). The raw data were deposited in the CNGB Sequence Archive (CNSA) of the China National GeneBank Database (CNGBdb, https://db.cngb.org/) under the accession number CNP0002115 (Guo et al. 2020). The adapter sequences were trimmed to remove low-quality sequences, and the high-quality clean reads were aligned to the reference genome of C. obscurus 7217 (CNGBdb accession number: CNP0001555) using hierarchical indexing for spliced alignment of transcripts (HISAT2) (Kim et al. 2015). Reads with no more than two mismatches from each sample were used to generate the transcripts using the StringTie software 1.3.1 (Pertea et al. 2016). The work ow of the HISAT-StringTie analysis is shown in Fig. S1B.

Quanti cation of the transcripts and functional annotation
Transcript expression was quanti ed as fragments per kilobase of transcript per million mapped reads (FPKM) using StringTie 1.3.1 (Pertea et al. 2016). Differential expression analysis of the transcripts between the three libraries was performed using the DESeq R package 1.10.1 (Trapnell et al. 2012). The resulting P-values were adjusted using the Benjamini-Hochberg method (multiple hypothesis test) for controlling the false discovery rate (FDR). The transcripts with an adjusted P-value of less than 0.01 and an absolute value of log 2 (fold change, FC) of more than 1 were designated as differentially expressed (Trapnell et al. 2012).

Co-expression network construction
The distinct hub genes of conidiation were established using the WGCNA (v1.69) package on BMKCloud platform (https://international.biocloud.net/) (Langmead and Salzberg 2012). Gene expression values were imported into WGCNA to construct co-expression modules using the automatic network construction function block-wise modules with the default settings. The expression levels of the differentially expressed genes (DEGs) were log-transformed using log 2 (FPKM + 1). Pearson's correlation coe cient was used to determine the co-expression relationship between each pair of genes. The WGCNA network was constructed with a soft thresholding power of β = 17, a minimum module size of 30 genes, the TOM-Type was unsigned, and the merge cut height was 0.25. The module-trait relationship was used to differentiate the hub genes among the conidiation stages.
Real-time quantitative polymerase chain reaction (RT-qPCR) RNA (1 µg) of each abovementioned and conidial samples was reverse-transcribed into cDNA using a PrimeScript RT reagent kit with gDNA eraser (TaKaRa). qPCR analysis of the cDNA samples was performed using SYBR Premix Ex Taq DNA polymerase II (TaKaRa), with the primers listed in Table S1. Experiments were performed on a qTOWER 2.2 real-time PCR thermal cycler (Analytik Jena, Germany) and the data was analyzed using qPCRsoft v1.1 (Analytik Jena, Germany). Transcript levels in the pre-C, C, post-C stages, and conidia were quanti ed using at least three independent replicates. The FC was normalized relative to the expression of the internal control gene encoding elongation factor 1-alpha (ef-1α), which was ampli ed using the primers EF1-F/R (Table S1), and calculated using the 2 −ΔΔCt method (Wang et al. 2018). The relative expression levels of the selected DEGs in the samples of the three stages were differentiated by one-factor analysis of variance (ANOVA) with the least signi cant difference (Fisher's LSD test) at a signi cance level of P ≤ 0.05. All analyses were conducted using an updated version of DPS software 7.05 (Tang and Feng 2007).

Global analysis of transcriptomes
The mycelia before, during, and after conidiation showed obvious morphological differences by changing from full and round hyphae (Fig. 1A), to spherical conidia formation in hyphal tips (Fig. 1B), and to shriveled hyphae (Fig. 1C), respectively. Approximately 90.7 Gb of clean reads were generated from the nine libraries of C. obscurus mycelia in the three stages, and 62.4-84.3% reads were mapped to the reference genome (Table S2). Transcriptome analysis identi ed 3,091 DEGs (1,879 upregulated and 1,212 downregulated) in C vs. pre-C, and 3,235 DEGs (1,954 upregulated and 1,281 downregulated) in C vs.

Function enrichment of the conidiation-related genes
The KEGG function enrichment analysis showed that the 3091 DEGs in C stage vs. pre-C stage were concentrated in amino sugar and nucleotide sugar metabolism, biosynthesis of amino acids, biosynthesis of antibiotics, carbon metabolism, and starch and sucrose metabolism ( Fig. 2A). The 3235 DEGs in C stage vs. post-C stage were concentrated in biosynthesis of amino acid, 2-oxocarboxylic acid metabolism, biosynthesis of antibiotic, and lysine biosynthesis (Fig. 2B).
The genes related to conidiation in C. obscurus were screened with the criteria of the genes upregulated in the C stage versus to either the pre-C or post-C stages. In total, 2,915 DEGs were putatively related to conidiation wherein 918 DEGs were upregulated in both C stage vs. pre-C stage and C stage vs. post-C stage (Fig. 3A). The KOG terms were concentrated in predicted general function, post-translational modi cation, protein turnover, signal transduction mechanisms, and amino acid transport and metabolism (Fig. 3B). The GO terms were concentrated in the following biological processes: singleorganism processes, response to stimuli, and signaling (Fig. 3C). Moreover, cellular components were enriched in membrane and membrane parts, whereas molecular functions demonstrated the most variations in catalytic activity, structural molecule activity, transporter activity, and electron carrier activity (Fig. 3C). The KEGG pathways were concentrated in the biosynthesis of amino acids and amino sugar and nucleotide sugar metabolism (Fig. 3D). A total of 375 DEGs putatively encoding transcription factors in 25 transcription factor families were predicted ( Fig. S3 and Table S3).

Co-expression network analysis
A WGCNA across all the nine samples elucidated the key conidiation-related genes. Highly interconnected genes were grouped into three different modules (black, brown, and turquoise) (Fig. 4). The black module containing 772 DEGs showed the strongest correlation with conidiation ( Fig. 4A) including 580 DEGs upregulated in both C stage vs. pre-C stage and C stage vs. post-C stage. Most of these DEGs were grouped into GO terms related to cell parts, catalytic activity, and metabolic processes ( Fig. 4B and Table   S4) and enriched in KEGG pathways involved in oxidative phosphorylation, ribosomes, and MAPK signaling ( Fig. 4C and Table S5). The KOG terms were concentrated in the following categories: posttranslational modi cation, protein turnover, chaperones, carbohydrate transport and metabolism, and cell wall/membrane/envelope biogenesis (Table S6).

Screening using the PHI database
The PHI database provides information on virulence factors produced in pathogen-host interactions. DEGs in the black module annotated by the PHI database showed that putative virulence factors were mainly associated with diverse hydrolases, including subtilisin-like serine proteases, trypsin-like serine proteases, lipases, and chitinases, which favor the invasion of pathogens into the host insect body through the cuticles (Table S7). A heat map using hierarchical clustering analysis showed that the nine samples localized into three clusters at each of the pre-C, C, and post-C stages wherein C stage DEGs had larger FPKM values than those in the other two stages (Fig. 5). Moreover, several DEGs were putatively related to secondary metabolite synthesis (Table S8), which may contribute to fungal virulence.

RT-qPCR assay
Four DEGs from the transcriptomic analysis were selected for further validation: heat shock transcription factor (EVM0008905), cytochrome P450 (EVM0004818), β-glucan synthesis-associated protein (EVM0003971), and trypsin-like serine proteinase (EVM0009373) using ef-1α as an internal control. The results from the three stages of conidia formation were consistent with the transcriptomic data (Fig. 6).

Discussion
The fungal conidia of entomopathogens play multiple roles in growth, dispersal, survival, and exploitation of nutritious resources in new habitats or hosts. Fungal conidiation is a key step in asexual reproduction, with morphological changes (tip growth) of conidiophores expanding in surface area and cell volume. This involves degradation, biosynthesis, and assembly (cross-linking between components) of conidial wall and membrane components such as glucan, chitin, and ergosterol. In this study, the genes encoding glycosidases, glucan synthases, chitin synthases, lipases, glucosidases, and those involved in mannoprotein and ergosterol synthesis and assembly were upregulated. The putative myosin chitin synthase gene (EVM0010183) may form bers of the cell wall (Weber et al. 2006), whereas glycosyl hydrolase (EVM0005673) is related to catalytic acidity, and 1,3-β-glucan synthase (EVM0007566) is related to cell wall remodeling and regulating cell wall thickness (Utsugi et al. 2002;Ishihara et al. 2007).
EVM0003314 encodes a putative transglycosidase homologous to the conidiation-speci c gene CRR1 and may be involved in spore wall assembly (Gómez-Esquer et al. 2004). EVM0001883 is putatively involved in intracellular tra cking, secretion, and vesicular transport and may mediate vacuole-vacuole fusion and cell membrane expansion (Tang et al. 2008). Ergosterol synthesis is essential for cellular growth and viability. This work revealed that the lipase ATG15 homolog (EVM0010219) was upregulated with previous work showing that it is involved in the regulation of intracellular sterol distribution and homeostasis (Ward et al. 2018).
C. obscurus conidia are multinucleated and contain dozens of nuclei (Barta and Cagáň 2006). The genes related to ribosome subunits, nuclear chaperones, and pentose phosphate pathways were upregulated (Tables S4 and S5), which may contribute to nuclear activities during conidiation. Substrate metabolism and energy metabolism was active during conidiation, and included the upregulation of genes related to the mitochondrial electron transport chain, which drives oxidative phosphorylation to produce adenosine triphosphate.
Many omics studies on the entomopathogenicity in Entomophthoromycotina revealed key fungal virulence factors including metalloproteases, lipases, subtilisin, and trypsin-like serine proteases (Grell et al. 2011;Małagocka et al. 2015;Wang et al. 2018). Most of them are secretory proteins, and high expression levels of the transcripts were found in fungal cultures on insect cuticle-inclusive medium, conidia, and mycotized cadavers (Małagocka et al. 2015;Keyhani et al. 2018;Wang et al. 2018). In this study, the genes encoding these secretory proteins were upregulated during conidiation without touching the host surfaces. This implied that conidiation in Entomophthorales is an auto-programmed event with the key virulence factors produced in advance to accelerate the infection. This may contribute to the higher virulence of Entomophthorales than that of Hypocreales against the same hosts (Pell et al. 2001).
Many DEGs putatively involved in secondary metabolite synthesis were upregulated during conidiation (Table S8). For example, EVM0003125 was upregulated 6-7-fold during conidiation; it is putatively related to plant-like isoquinoline alkaloid (fumisoquin) biosynthesis, which may lead to coloration and pathogenicity (Macheleidt et al. 2015;Baccile et al. 2016). EVM0001554 may mediate fumonisin biosynthesis which is a carcinogenic mycotoxin (Butchko et al. 2006). EVM0005600 may be involved in the biosynthesis of rubrofusarin (an orange polyketide pigment) (Rugbjerg et al. 2013). Con _423 may be involved in destruxin biosynthesis which is a cyclohexadepsipeptide involved in an immunomodulatory (Wang et al. 2012). EVM0003126 and EVM0008015 are putatively related to fungal trichothecene e ux pumps and may be involved in the synthesis of a atoxin precursor (Zhang et al. 2007;Bradshaw et al. 2009). EVM0005592 may be involved in the fusarin (a mycotoxin) biosynthesis (Niehaus et al. 2013). The tyrosinase homolog, EVM0007063 may be involved in the synthesis of the cyclic peptide ustiloxin B (Umemura et al. 2014). The potential synthesis of diverse toxic metabolites in C. obscurus implies that secondary metabolites may have an underestimated role in the virulence of conidia in Entomophthorales.
Conidia are aerial-born spores that have adapted to harsh conditions such as oxidative stress and UV radiation. Several DEGs related to environmental stress responses during conidiation were observed in this study. For example, proteins containing thiol-bearing cysteine residues that are sensitive to oxidation may contribute to oxidative stress resistance (Biteau et al. 2003). Copper/zinc superoxide dismutase is widespread throughout fungi and is proposed to promote catalysis through superoxide guidance (Gleason et al. 2014). The putative thioredoxin (EVM0006554) may be required for the repair of oxidatively damaged proteins and redox homeostasis (Trotter and Grant 2002). The copper/zinc superoxide dismutase homolog (EVM0004577) was expressed 3.72-fold during conidiation. Two upregulated DEGs (EVM0003905 and EVM0001948) homologous to ATP-dependent DNA helicase may be involved in DNA repair following UV-induced damage (Braga et al. 2015). EVM0001653 is putatively related to plasma membrane osmosensors that activate the high osmolarity glycerol (HOG) MAPK signaling pathway, which may be involved in the active ejection of conidia (Tatebayashi et al. 2006). Co_2537 is homologous to a two-component response regulator (SRR1) and may be required for stress adaptation, morphogenesis, and virulence (Desai et al. 2011).
The physiological activities involved in conidiation are complex and require transcriptional regulation. Long non-coding RNAs are involved in fungal virulence degradation in C. obscurus (Ye et al. 2021). In this study, many upregulated genes encoding diverse families of transcription factors were found in conidiation (Table S3 and Fig. S3). For example, EVM0003245 is homologous to the transcription factor Con7p in Magnaporthe grisea which regulates chitin synthesis in germinated spores and nuclear localization within spores (Odenbach et al. 2007). The C2H2 homeo-domain of zinc-nger transcription factors putatively involved in dimorphic transitions (Vallim et al. 2000) were found in many transcription factors during C. obscurus conidiation. A gene homologous to LaeA (EVM0004888) may play a pivotal role in inhibiting sexual development in response to light (Bayram et al. 2010). Moreover, EVM0002148 is related to the transcriptional regulation of genes responsible for the utilization of nonfermentable carbon sources (Turcotte et al. 2010). EVM0000205 is homologous to the transcription factor DEP6 that is involved in the biosynthesis of the histone deacetylase inhibitor depudecin, a virulence factor for Alternaria brassicicola (Wight et al. 2009). The putative heat shock transcription factor EVM0008905 is homologous to Prr1, which is involved in a His-Asp phosphorelay system in response to oxidative stress, cold temperature, and heavy metal toxicity (Ohmiya et al. 1999). EVM0002806 is homologous to stuA, and it is implicated in conidiophore development (Chung et al. 2015). It also implied that the diverse transcription factors in conidia favor the rapid response to surrounding environments.
In conclusion, transcriptomic pro ling revealed that conidiation of the aphid-obligate pathogen C. obscurus with diverse upregulated genes related to cell component biosynthesis, hydrolases, stress response, and secondary metabolite synthesis. This implied that the abundance of diverse transcripts is the basis for the multiple roles of conidia. These results provided a basis for further understanding the gene network for the infection cycle of entomopathogenic fungi and show the potential for conidiationbased formulations from entomophthoralean fungi as a natural pest-control.

Declarations
Authors' Contributions: XDZ and XZ conceived the idea, aims and designed the experiments. LZ, XW, TY and XZ did fungal culture, scanning electron microscopy, RNA extraction, and RT-qPCR., XZ did formal analysis. WY did data curation. XZ did writing-original draft preparation., XDZ did writing-review and editing. All authors have read through and agreed to the published version of the manuscript.
Funding: This research was funded by the grant from the National Natural Science Foundation of China (Grant Nos.: 31870637).
Availability of data and materials: Illumina sequence data have been submitted to CNGBdb (https://db.cngb.org/) under the accession number CNP0002155.
All data generated or analyzed during this study are included in this published article, accessions and its supplementary data le.
Adherence to national and international regulations: Not applicable.
Ethics approval and consent to participate: Not applicable. Consent for publication: Not applicable.
Competing Interest: The authors declare that they have no competing interest.  The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the differentially expressed genes (DEGs) in C stage vs. pre-C stage (A) and C stage vs. post-C stage (B). "Rich factor" refers to the ratio of the number of differentially expressed transcripts to the total number of transcripts for each pathway; the larger the rich factor, the higher the degree of enrichment. The q value indicated the P-value after the multiple hypothesis test corrections ranging from 0 to 1; the closer it is to 0, the more signi cant the enrichment considering a false discovery rate (FDR) ≤ 0.05 as the threshold.

Figure 3
Classi cation of upregulated genes in C stage vs. pre-C stage or C stage vs. post-C stage. (A) Number of upregulated genes putatively related to conidiation; (B) NCBI eukaryotic orthologous groups of proteins (KOG) functional annotation; (C) Gene Ontology (GO) functional annotation; (D) KEGG enrichment analysis. A smaller q value indicated a more signi cant enrichment, considering a FDR ≤ 0.01 as the threshold.

Figure 4
Weighted gene co-expression network analysis (WGCNA) identi cation of transcriptomes correlated with conidiation. (A) Hierarchical cluster tree of DEGs producing 3 gene co-expression modules; (B) The number of DEGs categorized by GO terms in the black module; (C) The number of DEGs categorized by KEGG pathways in the black module.

Figure 5
Heatmap of DEGs from the black modules of Figure 4 among the nine samples localized into three clusters at different stages (pre-C, C, and post-C). Genes with high and low FPKM values are orange and blue, respectively, whereas * refers to putative secretory proteins.