2.5.1. DNA Isolation

The isolation of genomic DNA was performed using the Genomic Mini AX Bacteria Spin kit (A&A Biotechnology, Gda ´nsk, Poland). All isolation steps were conducted according to the protocol provided by the manufacturer. The samples were taken after 168 h of the biodegradation process. The isolation efficiency was controlled by the fluorimetric method on the Qbit 3.0 device, using the Qubit™ dsDNA HS Assay Kit (ThermoFisher Scientific, Waltham, MA, USA). For each sample, three DNA extractions were carried out. Lastly, samples were mixed together, after a positive quantification.

## 2.5.2. NGS Sequencing

For microbial population taxonomy analysis, the V4 region of the 16SrRNA was amplified, based on the 515F-806R primers designed by Caporaso et al. (2012) [28]. The PCR details were optimized, and described in a previous publication [8].

Sequencing of the obtained amplicons was performed on the MiSeq platform (Illumina, CA, USA). The construction of amplicons and libraries normalization protocol is described in a previous work [26].

#### 2.5.3. Bioinformatic Analysis

The output sequencing data were analyzed using the CLC Genomics Workbench 8.5, with the CLC Microbial Genomics Module 1.2 software (Qiagen, Hilden, Germany). Readings were trimmed, demultiplexed, and paired ends were joined. Then, the chimeric readings were identified and removed. The data were clustered independently against two reference databases at 97% similarity of operational taxonomic units (OTU). The SILVA v119 database [29] was used for taxonomy annotation and biodiversity analysis, and the Green-Genes 13.5 database [30] was used for PICRUSt analysis. Selected biodiversity indices were determined: OTU number, Simpson's index, and phylogenetic diversity. Raw sequence data were deposited in the Sequence Read Archive (SRA) as project no. PRJNA831882.

In order to better characterize the analyzed microbiome, a linear discriminant analysis (LDA) effect size (LEfSe) [31] was performed. It allowed for the identification and selection of the most differentially abundant taxa. The microbial biomarkers that best describe the differences between the groups differing in the level of supplementation with nitrogen compounds were assessed. The following groups were determined: deficiency of nitrogen compounds: variants NP0 and NP1; and no deficiency (sufficient level: variants NP2 and NP3). LDA was coupled with effect size measurements: a non-parametric Kruskal–Wallis rank-sum test (*p* < 0.05), and the Wilcoxon rank-sum test (*p* < 0.05). A threshold of 3.5 for the logarithmic LDA score for discriminative attributes was established as the cut-off value.

#### 2.5.4. PICRUSt Analysis

In order to evaluate the bacterial functional composition in the analyzed samples, the bioinformatic tool PICRUSt v. 1.1.1 was used. This algorithm is widely used in environmental studies to analyze functional assessment of bacteria communities in different environments (soil, sediments, wastewater etc.) [32–38]. The predictions of functional composition of metagenomes were performed on the basis of sequencing data clustered against GreenGenes 13.5, and reference bacterial genomes deposited in the IMG database. The detailed algorithm's mode of action is described in a publication of its authors [39]. Based on the analysis, data on the predicted prevalence of 6009 KEGG Orthology IDs (KO IDs) participating in the degradation of polycyclic, aromatic hydrocarbons were obtained. It was assessed which OTU contributed to particular functions. Quality control of the PICRUSt prediction was performed according to the algorithm authors advices. The genome coverage was calculated for all samples using a weighted-Nearest Sequenced Taxon Index (NSTI) score [39].
