Identification of lncRNA and weighted gene coexpression network analysis of germinating Rhizopus delemar causing mucormycosis

ABSTRACT

Rhizopus delemar, an opportunistic fungal pathogen, causes a highly fatal disease, mucormycosis. Spore germination is a crucial mechanism for disease pathogenesis. Thus, exploring the molecular mechanisms of fungal germination would underpin our knowledge of such transformation and, in turn, help control mucormycosis. To gain insight into the developmental process particularly associated with cell wall modification and synthesis, weighted gene co-expression network analysis (WGCNA) was performed including both coding and non-coding transcripts identified in the current study, to find out the module of interest in the germination stages. The module-trait relationship identified a particular module to have a high correlation only at the resting phase and further analysis revealed the module to be enriched for protein phosphorylation, carbohydrate metabolic process, and cellular response to stimulus. Moreover, co-expression network analysis of highly connected nodes revealed cell wall modifying enzymes, especially those involved in mannosylation, chitin-glucan crosslinking, and polygalacturonase activities co-expressing and interacting with the novel lncRNAs among which some of them predicted to be endogenous target mimic (eTM) lncRNAs. Hence, the present study provides an insight into the onset of spore germination and the information on the novel non-coding transcripts with key cell wall–related enzymes as potential targets against mucormycosis.

KEYWORDS:

• Rhizopus delemar
• mucormycosis
• germination
• WGCNA
• Cazymes
• lncRNA
• eTM

1. Introduction

Mucormycosis (MCM), a serious and potentially fatal invasive fungal infection of humans, is caused by an opportunistic fungal pathogen called mucormycetes of the order Mucorales (Al Aboody and Mickymaray 2020). The disease has a high morbidity and fatality rate and usually establishes infection in immunocompromised individuals, particularly in those patients suffering from diabetes mellitus, bone marrow or organ transplantation, chemotherapy, neutropenia, corticosteroids treatment, elevated blood iron, haematological malignancy, and traumata (De Lira Mota et al. 2012). Currently, there are 27 different species from the order of Mucorales responsible for the mucormycosis disease and among them, the species Rhizopus delemar, previously known as Rhizopus oryzae, is the most common etiologic agent of mucormycosis and accounts for almost 70% of all cases (Skiada et al. 2020). The current treatment for mucormycosis is based on a combination of antifungal medications, with the dominant ones being- Amphotericin B (AmB), posaconazole, and isavuconazole and surgical removal of necrotic tissue, if necessary (Marty et al. 2016).

The fungal spores are the main causative agents of fungal infections, including mucormycosis, wherein germination is a crucial mechanism for developing infectious hyphae from the dormant spores. Although a healthy individual can withstand spore germination through phagocytic uptake, the mucoromycetes spores can survive within immune effector cells causing latent infections. In immunocompromised patients, failure of spore germination inhibition by phagocytes enables fungal growth and hyphal extension within tissues leading to angioinvasion, thrombosis, tissue necrosis, and ultimately, death (Manesh et al. 2016; Al Aboody and Mickymaray 2020). Therefore, studies to understand the critical developmental process of germination would expand our ways to explore the molecular mechanisms underpinning this transformation to control mucormycosis. Recently, a transcriptomic study conducted on Rhizopus delemar for a 24-hour time frame (from spore germination to hyphal growth) showed a shift in the gene expression pattern at different stages of germination where a high transcript-level was found to be associated with chitin-related processes at the dormant stage of germination (Sephton-Clark et al. 2018). Thus, the fungal cell wall is considered an essential aspect for the maintenance and survival of dormant spores and several studies have demonstrated the role of fungal cell wall organisation-related enzymes (chitin synthase, chitin deacetylase etc.) in maintaining its integrity (Ma et al. 2009; Xu et al. 2018). Moreover, the emerging role of long non-coding RNAs (lncRNAs) which are ≥ 200 nucleotides long sequences have been witnessed in response to cell wall remodelling, transcription regulation, nutrient response etc. in many fungal species (Hovhannisyan and Gabaldón 2021). These lncRNAs based on their mechanism of interference are categorised as cis or trans acting (Camblong et al. 2007, 2009) and recently their elaborated functional role in various categories have been addressed in human diseases (Chen et al. 2021). But in most of the fungi, the functional lncRNAs identified are cis acting especially interfering with the transcription of closely situated genes (Li et al. 2021), such as nc-tgp1, prt, and prt2 (Shuman 2020). Still, knowledge of lncRNA on the fungal species including Rhizopus delemar, is very limited and perhaps little is known with regard to its role on R. delemar virulence and mucormycosis disease.

RNA-Seq is a widely used high-throughput gene expression profiling technology to provide an efficient and descriptive understanding on the molecular basis of a given sample over time through differential gene expression. However, to analyse these large datasets for molecular understanding of biological processes at a system level, network approaches prove to be an attractive approach (Barabási and Oltvai 2004). As genes with similar expression patterns can be associated with their functional roles, identifying these correlated genes helps to shed light on their possible roles in biological processes. Weighted gene co-expression network analysis (WGCNA) also known as differential co-expression network analysis is a systems biology approach that clusters genes with similar co-expression profiles into several characteristic modules. This approach by applying the scale-free topology criteria and having co-expression similarity raised by β > 1(weight) prevents information loss compared to the conventional unweighted co-expression analysis approach (Zhang and Horvath 2005; Langfelder and Horvath 2008). Thus, based on the previous transcriptomic profiling of R. delemar germinating spores, the present study aims to understand the germination process with respect to cell wall modifying enzymes in a more holistic way using the weighted gene co-expression approach including the identification and co-expression of novel lncRNAs that may have a role in controlling the expression of the important cell wall-related enzymes.

2. Materials and methods

2.1. Dataset retrieval

The transcriptomics data of a germinating Rhizopus delemar strain RA 99 880 was downloaded from National Centre for Biotechnology Information (NCBI) – Gene Expression Omnibus (GEO) database with study accession PRJNA472797. The dataset included 30 paired-end samples sequenced (Supplementary Table S1) on the Illumina HiSeq2000 platform at different time points of germination. These were named as resting phase for 0 hour, isotropic phase for 1–6 hours and hyphal growth phase for 12, 16, and 24 hours for this study, with three biological replicates at each time point (Sephton-Clark et al. 2018).

2.2. Data pre-processing and identification of unannotated transcripts

All the RNA-seq raw data were uploaded on Galaxy suite (https://usegalaxy.org/) and an initial quality check was done using FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The short reads were then mapped to R. delemar reference genome (RA 99–880 assembly RO3) using HISAT2 (Kim et al. 2015) for generating the binary alignment matrix (BAM) files following which MultiQC was performed. Reads quantification and novel transcript identification were performed using StringTie v.2.1.1 (Pertea et al. 2016). For this, individual transcript assemblies were generated and combined using StringTie’s merge option. Further, the merged output in GTF format was compared with the reference gene annotation file using gffcompare (Pertea and Pertea 2020) to classify transcripts into different classes. Finally, the unannotated transcripts with class codes “u” (intergenic), “x” (anti-sense), “i” (intronic), “o” (generic exon overlap with reference transcript), and “e” (single exon transfrag overlapping reference exon) were extracted.

2.3. Prediction of novel lncRNA

For the identification of lncRNAs, the following steps were performed sequentially. Firstly, the transcripts with the number of exons ≥ 2 and sequence length ≥ 200 were filtered and sequences were extracted using the GFFread tool by taking the reference genome file and GTF file of the filtered transcripts. Next, the identified transcripts, if they had any coding potential, were screened using PLEK (Li et al. 2014) and FEELnc (Wucher et al. 2017) and the putative lncRNAs were filtered out from the intersection of non-coding transcripts identified from these two methods. Subsequently, the predicted lncRNA transcripts with at least one significant hit (e-value: 1e⁻⁵) against Pfam (Mistry et al. 2021), Swiss-Prot and TrEMBL (The UniProt Consortium 2021) using blastx in NCBI standalone tool and Blast2GO (Götz et al. 2008) were further discarded. Finally, tRNA was scanned through tRNAscan-SE (http://lowelab.ucsc.edu/tRNAscan-SE/) followed by an rRNA scan by Rfam blast and the remaining transcripts after filtrations were considered as putative lncRNAs.

2.4. Weighted gene co-expression network analysis

To construct a co-expression network and identify significant modules at different stages of R. delemar germination, the “WGCNA” package on R was used (Langfelder and Horvath 2008). The read counts were generated using StringTie eB DEseq2 count generation pipeline, which was further normalised using the DESeq2 tool (Love et al. 2014) variance stabilising transformation (VST) and visualised through principal component analysis (PCA) plot. To proceed with weighted adjacency matrix creation following a scale-free network, the suitable soft threshold power (β) and its corresponding mean connectivity were calculated using the “pickSoftThreshold” function in WGCNA (Zhang and Horvath 2005). Subsequently, the adjacencies calculated were transformed into Topology Overlapping Matrix (TOM) to determine the network connectivity and the value of dissTOM (1-TOM) was used for module detection using unsupervised hierarchical clustering corresponding to clusters (modules) having similar gene expression with minimum module size = 30, followed by DynamicTree cut algorithm where modules with similar expression patterns were merged and visualised. Next, the module-trait association was done with traits defined as resting phase, isotropic phase and hyphal growth phase and a heat map of the module-trait relationship using Pearson’s correlation coefficient and an expression heatmap for selected module was generated. Finally, the module with the highest correlation was identified as a key module in this study and module membership (MM) vs. gene significance (GS) of the corresponding module was further analysed.

2.5. Annotation and gene enrichment analysis of significant module

The transcript accessions from the identified key module were further filtered out in the GS/MM correlation. These final transcripts were taken for gene ontology (GO) annotation, carbohydrate active enzyme (CAZyme) annotation and gene enrichment analysis. For GO annotation, the Blast2GO suite (Götz et al. 2008) was used with an inbuilt InterproScan tool by selecting default parameters. The GO annotation of the selected module was visualised using the WEGO tool (https://wego.genomics.cn/) and the enrichment analysis was done on AgriGO V2.0 (Tian et al. 2017) using single enrichment analysis (SEA) and further assessed through Revigo (Supek et al. 2011). For CAZyme annotation, dbCAN2 tool (http://cys.bios.niu.edu/dbCAN2) was used with default parameters and only those CAZymes were retained which were predicted through more than two integrated tools of dbCAN2.

2.6. Network construction and lncRNA-miRNA-mRNA interaction

From the selected module, nodes were extracted with an edge weight cut-off of 0.5 using the function “exportnetworktoCytoscape” in the WGCNA package. The network was visualised and annotated for lncRNA and mRNA on Cytoscape (Shannon et al. 2003) and sub-clustering analysis was done using the Cytoscape plugin MCODE (Bader and Hogue 2003) with default parameters. Further, among the accessions in the network, mRNA-lncRNA interaction was identified using LncTar (Li et al. 2015) which finds the interaction using free energy minimisation with normalised binding free energy (ndG) cut-off =< −0.01. Further to find out the lncRNAs in the network that could act as an endogenous target mimic (eTM) of miRNA, psMimic tool was utilised (Wu et al. 2013). For this, all the mature miRNA fasta format files were downloaded from miRBase database (https://mirbase.org/) and clustered using CD-HIT-EST tool (Fu et al. 2012) with default parameters before predicting the eTM lncRNA. The miRNAs predicted to have mimicry with lncRNA, were used to identify the target mRNA within the network using psRobot tool (Wu et al. 2012) with moderate settings.

3. Results

3.1. Identification of novel lncRNAs

Based on the pipeline (Figure 1), the fastq files after mapping with the reference genome, resulted in an average of 82% read alignment (Supplementary Figure S1). The total number of assembled transcripts including both existing and novel with class codes from the merged GTF file resulted in 39,533 transcripts [Supplementary Data D1, 10.6084/m9.figshare.24099387); this file can be used for extracting fasta sequences from the genome sequence of R. delemar, RA 99–880 assembly RO3, using gffcompare tool. Among them 17,503 were completely matched (class code “=”) and from the remaining transcripts, 6,787 transcripts belonged to class codes (“x”, “i”, “o”, “e”, “u”). The class codes distribution of the merged transcripts is provided in Supplementary Data D2]. These 6,787 transcripts after filtration based on the number of exons ≥ 2 and length of transcript ≥ 200 bp, were further reduced to 6,571 transcripts. Among these, 4,498 transcripts belonged to class code “u” (intergenic), 1,458 belonged to “x” (overlap with a known gene on the opposite strand or anti-sense), 609 belonged to “o” (overlap with known exon), six belonged to “i” (intronic), and zero transcripts in “e” (single exon transfrag with partial intron covering). The length distribution and exon number were also calculated for the identified lncRNAs (Figure 2a), where the sequence length ranged from 201 bp to 3,538 bp, with 72.88% in the range of 200–1,000 bp and 27.12% were found in more than 1,000 bp. The exon number distribution revealed that nearly more than half of the lncRNAs consisted of two exons, and only 11 lncRNAs had exon number > 10.

Figure 1. The overall workflow of lncRNA identification and weighted gene co-expression network analysis.

Figure 2. (a) Distribution of non-coding transcripts after cut-off ≥2 exons, ≥ 200bp; (i) Length distribution of putative lncRnas; (ii) Exon number distribution. (b) Venn diagram of (i) Potential lncRnas predicted from two tools (PLEK and FEElnc) for identification of commonly annotated lncRnas; (ii) 4,800 non-coding transcripts for potential coding hits identified through blastx from the three databases, viz. Pfam, Swiss-Prot and TrEMBL. (c) Distribution of final non-coding transcripts after filtration steps; (i) Length distribution of majorly represented class codes, i.e. intergenic (u) and anti-sense (x); (ii) Exon number distribution of intergenic (u) and anti-sense (x) transcripts.

Further, transcripts designated as lncRNAs were screened for coding potential and removed by comparing the commonly predicted non-coding transcript accessions among PLEK (5166) and FEELNc (5826), resulting in 4,800 non-coding transcripts. These transcripts were further subjected to blastx against three customised databases using nucleotide sequences from Swiss-Prot, Pfam and TrEMBL combined which resulted in 2,127 accessions with blast hits (Figure 2b) and 2,673 unique accessions without blast hits. In the next filtration step for tRNA and rRNA scan through tRNAscan-SE and Rfam blast, respectively, nine rRNAs were detected and no tRNA was found. These filtered lncRNAs (2,664) were further subjected to blastx on Blast2GO that ended up with 95 accessions that could have coding potential recognised by PROSITE and mobidb-lite (identify intrinsically disordered regions) (Figure 2b). Finally, a total number of 2,569 potential lncRNAs were obtained that were majorly classified as intergenic (u), and antisense (x) lncRNAs, i.e. 2,346 and 218, respectively, and five lncRNAs classified as intronic (i) based on their genomic locations. Among these final lncRNAs, a similar trend was observed in the sequence length distribution, where nearly 50% of the sequences were in the range of 200–500 bp in each class with a maximum length of 2,025 bp in “u” and 1,722 bp in “x” class codes (Figure 2c).

3.2. Key module identification by WGCNA

The read counts that were generated for each time point from DESeq2 using the merged GTF file and the individual BAM files through StringTie -e and prepDE.py python code revealed 41% of variance through PC1 and 19% by PC2 through the PCA plot of VST normalised read counts of the overall samples (Figure 3a). The PCA in this study demonstrated no batch effects among the samples in which all three phases of germination, i.e. resting phase, isotropic phase, and hyphal growth phase were grouped in distinct clusters. This was coherent with the TMM-based PCA of the original author’s work by Sephton-Clark et al. (2018). Additionally, the network topology was screened to construct the gene co-expression network by choosing a suitable soft threshold underlying a scale-free network. Therefore, the soft threshold power value of 18 (β = 18) having a correlation coefficient (R²) >0.8 was chosen because of its corresponding mean connectivity of 101 (Supplementary Table S2, Figure 3b). As for the given soft power based on R², the selection was preferred for the lower mean connectivity. Next, hierarchical tree clustering was performed based on the adjacency matrix and its subsequent topology overlapping matrix between the genes. Thus, gene modules with similar expressions were assigned with distinct colours and genes unassigned to any module were represented by grey module (Figure 3c). Further, modules with high expression similarities were merged using module eigengene (ME) by keeping cut-off (MEDissThres = 0.25) that generated a total of 40 co-expressed modules based on the gene expression similarities (Figure 3d, Supplementary Table S3) which were then associated with the sample traits, i.e. the different phases of germination. Thus, a heat map representing module-trait association is shown in Figure 4a, where blue to red indicates module trait correlation from negative to positive. The module-trait association analysis revealed a particular (purple) module to be significantly correlated (Pearson correlation = 0.95) only in the resting phase having 4,027 accessions. Moreover, other positively correlated (Pearson correlation > 0.7) modules with specific traits were also found, such as maroon module having 64 accessions in the first hour isotropic phase. Similarly, modules lavenderblush3 and darktorquiose with 52 and 593 accessions in each were correlated with the hyphal growth phase at 24 hours. Therefore, the highest positive correlation (0.95) between the purple module and resting phase, i.e. onset of germination was chosen to be our module of interest. Hence, an expression heat map for the purple module corresponding to all the phases was also generated to visualise the expression pattern among the three major phases of germination i.e. resting, isotropic, and hyphal growth, which revealed a distinct expression pattern in the resting phase (Figure 4b). This was further assessed through correlation analysis between gene significance (GS), which is a correlation between a trait and gene expression profiles, and module membership (MM), which is a correlation of a gene expression with module eigengene, i.e. (GS vs MM) of the purple module in each of the development stages that showed a high correlation value of 0.8 and P-value <1e⁻²⁰⁰ in the resting stage (Figure 4c). Finally, the list of correlated genes present in the purple module for the resting stage was further filtered out with p-value cut-off 0.05 individually for both MM and GS in module annotation This filtration resulted in a total of 3,787 accessions in purple module, from which 425 accessions belonged to lncRNAs and within the remaining 3,362 accessions of coding transcripts, consisted 1,081 existing and 2,281 novel transcripts.

Figure 3. (a) 2D PCA plot of normalised count samples explained by 41% variance on the x-axis and 19% variance on the y-axis with clear distinction in the sample clustering. (b) Scale-free and mean connectivity with soft power values on the x-axis and correlation coefficient of model fit and corresponding mean connectivity represented on y-axis. (c) Gene dendrogram with assigned module colours. (d) Merged modules after cut-off: 0.25 for similarly expressed gene modules (grey modules consist of genes with no similarities in the expression).

Figure 4. (a) Correlation heat map of module-trait association for 40 modules, ranging from red to blue, indicating positive to negative correlation with p-value mentioned in the bracket for the defined traits in each module. (b) Expression heat map of the purple module for each trait. (c) Correlation of module membership (MM) of purple module on x-axis and gene significance (GS) of each phase on y-axis i.e. resting (i) isotropic (ii–vii) and hyphal growth (viii–x). (R: resting, I: isotropic, H: hyphal).

3.3. Key module annotation and enrichment

Within the purple module for the resting stage, the GO annotation resulted in 1,607 accessions associated with GO information, from which 1,010 were novel transcripts. Additionally, 842 accessions had domain information and classification from various databases (Supplementary Data D3), of which 527 were reported as novel transcripts. The GO ids of the accessions after visualisation (Figure 5a) revealed majority of the accessions belonged to the membrane, cell and cell parts performing catalytic and binding activities involved in metabolic and cellular processes followed by localisation processes. These accessions were further taken for their statistical significance, which through hypergeometric test and FDR cut-off ≤0.05 revealed the accessions in the purple module were enriched for phosphorus and phosphate-containing metabolic process (GO:0006793, GO:0006796), transmembrane transport (GO:0055085), regulation of various metabolic processes, phosphorylation (GO:0016310), organic cyclic compound biosynthetic process (GO:1901362) and cellular response to stimulus (GO:0051716), further among these processes it was revealed that protein phosphorylation, carbohydrate derivate metabolic process (GO:1901135) and cellular response to stimulus to be more specific processes. Moreover, the functions for which the module was enriched were mainly binding functions (GO:0005488), including heterocyclic compound (GO:1901363), organic cyclic compound (GO:0097159), protein binding (GO:0005515) and DNA binding (GO:0003677). Other enriched functions were catalytic (GO:0003824), oxidoreductase (GO:0016491), transferase (GO:0016772, GO:0016740), kinase (GO:0016301), DNA binding transcription factors (GO:0003700), hydrolase (phosphoric ester) (GO:0016788), transferase (GO:0016740), hydrolase (ester) (GO:0016788), transmembrane transport (GO:0022857), and transport activities (GO:0005215) (Supplementary Table S4, S5) and the more specific functions were hydrolase, DNA binding transcription factor, protein binding, kinase, and transferase (phosphorus-containing groups). Further, by taking the nucleotide sequences of the transcript accessions in the module for CAZymes annotation, a total of 63 CAZymes were annotated after filtering the prediction with ≥ 2 tools. Among the 63 CAZymes, seven belonged to families of auxiliary activities (AA), 33 to glycosyl hydrolases (GH), 19 to glycosyl transferases (GT) and one each in carbohydrate-binding modules (CBM), carbohydrate esterase (CE) and two in polysaccharide lyases (PL) (Figure 5b, Supplementary Table S6). Apart from one family each in CE and CBM modules, most families represented in these groups were AA1 (laccase), GT2 (chitin synthase), GT15 (mannosyl transferase), GH28 (polygalacturonase), GH16 (chitin β-1,6-glucanosyltransferase), GH18 (chitinase) and PL38 (glucuronan lyase) disclosing their role in cell wall synthesis process.

Figure 5. (a) Gene ontology (GO) distribution of the annotated genes of the purple module with the overrepresented GOs labelled on the x-axis and relative percentage/numbers for each classification on the y-axis. (b) Number of CAZymes in each modules AA, CBM, CE, GH, GT, and PL. (c) Distribution of CAZyme families in each module. AA: auxiliary activities, CBM: carbohydrate binding modules, CE: carbohydrate esterases, GH: glycosyl hydrolases, GT: glycosyl transferases, PL: polysaccharide lyases.

3.4. Network analysis and LncRNA target prediction

The purple network exported from the WGCNA module with an edge weight cut-off of 0.5 resulted in 264 nodes with 11,817 edges. These nodes comprised 102 accessions with GO and InterPro IDs, 65 predicted lncRNAs and 97 accessions with no known conserved domains, with most of them generally coding for short amino acid sequences. Among the 102 annotated nodes, six were coding for CAZymes, viz. RO3G_03340 (GH16_19, chitin β-1,6-glucanosyl transferase), MSTRG.5422.2 (GH28, polygalacturonase), RO3G_09825 (GH46, chitosanase), RO3G_06658 (GH47, α-mannosidase), RO3G_04895 and RO3G_01689 both coding for GT15 (glycolipid 2-α-mannosyltransferase). The protein-coding accessions having GO information after hypergeometric test and Hochberg FDR cut-off ≤0.05 were enriched for the functions of ion binding activity (majorly zinc ion transport, GO:0043167), substrate-specific transporter (majorly solute: protein antiporter, GO:0022892) activity and hydrolase activity (majorly on carbohydrate substrates, GO:0016798) in the network.

Further, MCODE clustering analysis resulted in two sub-clusters, with sub-cluster 1 (cut-off: 100.013) having 157 nodes with 7801 edges and sub-cluster 2 (cut-off: 3.846) having 14 nodes with 25 edges. The sub-cluster 1 consisted of four CAZymes, 41 lncRNAs and 112 protein-coding accessions, whereas sub-cluster 2 consisted of one CAZyme, three lncRNAs and ten protein-coding accessions. The CAZymes present in the subcluster 1 were GT15 (glycolipid 2-α-mannosyltransferase), GH28 (polygalacturonase) coded by a novel transcript (MSTRG.5422.2), GH46 (chitosanase), and GH47 (α-mannosidase). Moreover, the nodes among the sub-cluster 1 with > 200 degrees were GT15 (glycolipid 2-α-mannosyltransferase) with a degree of 226 and 11 lncRNAs (MSTRG.12565.1, MSTRG.10396.1, MSTRG.2073.5, MSTRG.4963.2, MSTRG.11048.1, MSTRG.8879.1, MSTRG.4438.1, MSTRG.9949.1, MSTRG.1705.2, MSTRG.776.1, MSTRG.4508.1) with degrees ranging from 233 to 216, including 19 other unannotated coding transcript accessions. In sub-cluster 2, the only CAZyme present was GH16_19 (chitin β-1–6-glucanosyl transferase) and three lncRNAs (MSTRG.1226.2, MSTRG.2874.1, MSTRG.4519.2) of which lncRNA accession MSTRG.2874.1 was having co-expression with most of the nodes in this sub-cluster (Figure 6). Further, among the accessions present in the network, lncRNA-mRNA interaction study revealed 62 lncRNAs to be interacting with 131 coding mRNAs with normalised free energy ranging from −0.1001 to −3.29 (Supplementary Table S7), out of which 5 lncRNAs were acting as antisense lncRNA having interaction with 30 coding accessions (17.96%) in the network. Among the CAZymes, except for GH46 (chitosanase, RO3G_09825) all the CAZymes were interacting with lncRNAs i.e. GH16 (chitin β-1–6-glucanosyl transferase, RO3G_03340) interacted with lncRNAs (MSTRG.10396.1 and MSTRG.13812.1) which were one of the top degree distributed nodes in subcluster 1, GH28 (polygalacturonase, MSTRG.5422.2) with MSTRG.7652.1, GH47 (α-mannosidase, RO3G_06558) with MSTRG.220.5, GT15 (α-mannosyltransferase) coded by two accessions RO3G_01689 and RO3G_04985 interacted with MSTRG.10396.1 and MSTRG.13812.1, respectively (Table 1). Moreover, among the top 11 lncRNAs in the subcluster 1 co-expression network, except for lncRNA (MSTRG.2073.5), all were having mRNA interaction, whereas in subcluster 2, all the three lncRNAs were interacting with the coding mRNAs in the main network. This suggested the highly co-expressed network exported with an edge weight cut-off of 0.5 also had important lncRNA interactions within the network implicating their role as positive trans-acting regulators.

Figure 6. (a) Overall network representation of the 264 nodes exported from WGCNA with edge weight cut off 0.5, in degree attributed grid layout, where S1 and S2 are the network subclusters and the grid nodes on the right side involved in lncRNA-mRNA interaction with normalised binding energy (ndG) =< −0.1. (b) Highly interconnected sub-cluster 1 (MCODE score: 100.013) from the main network with the labelled lncRNAs and CAZyme accession nodes that were having >200 degrees. (c) Subcluster 2 (MCODE score: 3.846) showing one CAZyme (GH16_19) and three lncRNAs. GH: glycosyl hydrolase, S1: Subcluster 1, S2: Subcluster 2.

Table 1. Interaction of lncRNA and CAZymes within the network exported from purple module.

CAZymes	Name	Accession	Interacting lncRNA	dG	nDG
GH16_19	Chitin β-1,6-glucanosyl transferase	RO3G_03340	MSTRG.10396.1	−13.07	−0.2042
			MSTRG.13812.1	−12.35	−0.1187
GH28	Polygalacturonase	MSTRG.5422.2	MSTRG.7652.1	−9.43	−0.1014
GH47	α-mannosidase	RO3G_06558	MSTRG.220.5	−7.53	−0.1345
GT15	α-mannosyltransferase	RO3G_01689	MSTRG.10396.1	−10.2	−0.1437
		RO3G_04985	MSTRG.13812.1	−11.78	−0.1042

dG: Free energy, nDG: Normalised free energy i.e. dG/min (length lncRNA; length mRNA).

Additionally, in the network, 32 lncRNAs were identified as eTM of 25 miRNA (Supplementary Table S8) and among these miRNAs, 6 miRNAs were identified to have targets in the network i.e, ath-miR407 targeting MSTRG.1594.1 (an uncharacterised novel accession), bta-miR-2325c targeting MSTRG.11071.2 coding for sodium-proton anti-porter activity, MSTRG.6829.2 and MSTRG.6829.3 both had no known conserve domain revealed through protein-blast, hsa-miR-302f targeting RO3G_13192 coding for a hypothetical protein, has-miR-4306 targeting RO3G_14505 coding for smc superfamily activity (COG1196, cl34174) that involves in chromosome segregation ATPase activity controlling cell division and chromosome partitioning, hsa-miR-4306 targeting MSTRG.13956.2 which codes for major intrinsic protein (MIP) superfamily (CDD id: cl00200) involved in channel activity and stu-miR-172c-5p targeting MSTRG.2016.2 which through blastx revealed aromatic acid exporter family member 2 superfamily, specific to fungi and yet to be functionally characterised (Supplementary Table S9).

4. Discussion

Rhizopus delemar, which created havoc during the COVID-19 pandemic, is an opportunistic fatal fungal pathogen, but the pathogenicity mechanism is still in infancy. For successful pathogenesis, the onset of germination, i.e. activation of the resting spore stage, is a crucial phase which constitutes a series of functional changes, including transcriptional activation, process of cell wall synthesis and remodelling. In the current study, the annotation and enrichment analysis of the co-expression module in the dormant stage also revealed the other activities, apart from the activities such as oxidoreductase, hydrolase and transferase activities, transfer of phosphorus-containing group as reported by Sephton-Clark et al. (2018). These activities were mainly protein binding activities belonging to various processes such as signal transduction, intracellular protein transport, ubiquitin-related, and protein prenylation. Moreover, DNA-binding transcription factor activities involved in DNA integration and DNA templated transcription regulation were also enriched hence unfolding the mechanism of onset of active division in spore cells.

The fungal cell wall is the first structural component that comes into contact with the host during an infection or spore germination and acts as a key factor in host-pathogen interaction. The spore wall is composed of melanin, N-acetylglucosamine, chitosan, mannan, and glucans (Bartnicki-Garcia and Reyes 1964), and through cell wall remodelling grows into hyphae, where the enzymes involved in this process can be the key candidates for drug target. The temporal expression of these cell wall-modifying enzymes during the germination process in Mucorales still has a knowledge gap (Lecointe et al. 2019). The fungal cell wall synthesis process, including chitin biosynthesis, is an imperative process and the role of chitin synthase (CHS) belonging to the CAZyme family glycosyltransferase 2 (GT2) has been well established as a crux for chitin biosynthesis process (Coutinho et al. 2003). As plants and animals (humans) lack chitin, CHS serves as an important target for antifungal agents. Moreover, CHSs are directly associated with hyphal growth and many studies demonstrating the association of reduced pathogenic virulence with inactivated CHS genes have also been reported (Nix et al. 2009; Chaudhary et al. 2013; Li et al. 2019). The identification of the GT2 (chitin synthase) family in abundance at the initial stage of germination (dormant phase) in the current study further supports its significance. It is important to note that fungal cell wall plasticity and maintenance requires continuous remodelling of the cell wall that involves the processes of synthesis, crosslinking and degradation. The identification of a high copy number of the GH18 family (chitinase) co-expressed with chitin synthase at the dormant stage of fungal germination witnessed in this study also corroborates with the previous work showing a high expression level of chitinases in ungerminated R. delemar (Sephton-Clark et al. 2018) and Aspergillus niger (van Leeuwen et al. 2013) spores. Moreover, polyphenol oxidation contributes to the pathogen virulence (Zhu and Williamson 2004) and the identification of CAZyme AA1 (laccase) with three copy numbers also has a significant role due to their involvement in polyphenol oxidation, represented in high copy number within the auxiliary activities (AA) families in the resting phase. Hence, the role of AA1 in mucor spore germination could be another interesting domain apart from well-studied enzymes such as chitin synthase and chitinase.

The fungal lncRNAs have diverse functions regulating various biological processes, including fungal development (Donaldson and Saville 2012). The lncRNAs are known to modulate vegetative growth, phosphate regulation and cell-to-cell adhesion in yeast (Li et al. 2021) and are also involved in the regulation of cellulase genes (HAX1) in Trichoderma reesei (Till et al. 2018). Still, the regulatory role of lncRNA in the pathogenic fungi is lacking (Dhingra 2020), hence the present study aimed to identify lncRNA in Rhizopus delemar and their regulatory roles. Therefore, to find out the role of lncRNA in regulating the cell wall–related enzymes and other important enzymes through co-expression and lncRNA-mRNA interaction, the study revealed 11 highly interconnected novel lncRNAs that were co-expressed with accessions coding for cell wall degrading and mannose sugar utilising CAZymes present in the network sub-cluster as well as their interaction with other coding mRNAs in the network suggested their important role in positive regulation. Among the CAZymes accessions in the network, apart from their co-expression, all of them were found to be interacting with lncRNAs except for GH46 (chitosanase, RO3G_09825) acting on chitosan, a deacetylated product of chitin. Interestingly, a CAZyme GH28 (polygalacturonase) encoded by a novel transcript MSTRG.5422.2, was also co-expressed with all the 11 highly interconnected lncRNAs and found to be interacting with one lncRNA (MSTRG.7652.1) in the network. This family is an important cell wall degrading enzyme (CWDE) which is involved in pectin metabolism and crucial for pathogenesis when infecting a plant host by utilising pectin material for its growth and pathogenesis (Kubicek et al. 2014); however, its role in fungal cell wall modification or spore germination is still elusive which gives a scope to explore more on this family.

Further, two CAZymes in network subcluster 1, GH47 (α -mannosidases) mainly engaged in the remodelling of polysaccharide portion of the fungal cell wall that ensures its robustness and maintenance (Munro et al. 2005) and GT15 (glycolipid 2-α-mannosyltransferase), involved in mannose utilisation were co-expressed and had interactions with lncRNA with latter (RO3G_04985) having a high degree distribution. Apparently, mannose is an important constituent of the fungal cell wall and the enzyme glycolipid 2-α-mannosyltransferase (GT15) involved in the synthesis of polysaccharides and other glycoconjugates constituting mannose could be a prospective drug target (Wagener et al. 2008). Moreover, studies demonstrating the disruption or deletion of the mannosyltranferase gene affecting fungal adhesion, hyphal development and eventually attenuation of fungal virulence have also been reported and suggest it as a potent antifungal target (Fabre et al. 2014; Díaz-Jiménez 2017). Another key CAZyme, apart from GT15, GH28, and GH47 in the highly interconnected sub-networks was GH16 which was found to be co-expressing with the lncRNA (MSTRG.2874.1) and interacted with two lncRNAs (MSTRG.10396.1 and MSTRG.13812.1). This CAZyme GH16 subfamily 9 (chitin β-1–6-glucanosyl transferase), is an enzyme required for strengthening of the fungal cell wall during cell wall biogenesis through chitin-glucan (β-1–6) crosslinking (Arroyo et al. 2016) which is crucial for cell wall integrity as it was evident that the Rhizopus species lacks β-1,3-glucans (Ostrosky-Zeichner et al. 2005). Therefore, this enzyme could also be further utilised as a potential candidate for drug target given its role in fungal cell wall integrity. Moreover, the lncRNAs identified in the network acting as eTM of miRNA for negative regulation revealed further that these miRNAs were involved in controlling transport activities especially solute transport and cell division which are crucial for cell growth and viability.

Overall, the initial events of spore germination apart from revealing the importance of signal transduction process and post translational modification for intracellular transport, also revealed the well-established role of cell wall–related enzymes such as chitin synthase and chitinase. Moreover, the identification of novel lncRNAs and several other cell wall–related enzymes, such as polygalacturonase, mannosidase, mannosyltransferase and transglycosylase and their interaction with novel lncRNA suggested their role to be crucial in the process of germination.

5. Conclusion

The mucormycete spores exist abundantly in the environment and are capable of infecting immune-compromised patients. The present study, apart from annotating novel coding and non-coding transcripts, also revealed some of the important cell wall-related enzymes crucial for the onset of germination. Further, the co-expression network analysis identified the highly clustered nodes representing CAZymes and lncRNAs. Among these CAZymes some of them were also represented high copy numbers in the significant module of the resting phase. Most of the highly interconnected lncRNAs identified were co-expressed with candidate CAZymes involved in mannosylation, chitin-glucan crosslinking, polygalacturonase and chitosanase activities as well as had trans-acting positive regulation. Therefore, this study provides a knowledge base and a future scope to target and characterise these hub genes especially CAZymes and the identified lncRNAs. Further, elucidating the role of lncRNAs in the expression of these CAZymes and other crucial enzymes present in the co-expression network would help pave the way to combat mucormycosis.

Acknowledgments

Barsha Kalita is supported by Indian Council of Medical Research (ICMR), New Delhi, India, in the form of ICMR-Senior Research Fellow (2021-10524). The authors also acknowledge the library facilities extended by Pondicherry University.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/21501203.2023.2265414