Microplastics are a ubiquitous presence in the world’s aquatic environments and their threat to aquatic biota is poorly understood, especially in freshwater ecosystems. In the environment, microbial biofilms can form on the surface of microplastics, and these plastics have the potential to adsorb harmful toxins. Because lab-based studies on microplastics are often conducted with clean polymers, in ecologically unrealistic conditions and concentrations, the impact of these weathered microplastics on aquatic organisms in ecologically realistic conditions is still unclear. To help address the need for ecologically relevant microplastic exposure data, we incubated 500 μm polyethylene microplastic beads in Muskegon Lake, Michigan, USA and used them to conduct a 28-day ingestion study with male and female fathead minnows (Pimephales promelas). We examined the effects of microplastic ingestion on the fish gut microbial community along with hepatic gene expression and health parameters. We found that microplastic ingestion had statistically significant impacts on growth in male fathead minnows. Microplastic treatment did not significantly alter the beta diversity of the gut microbial community for either males or females, but there were clear differences between sexes and over time, indicating that these factors may outweigh the impacts of microplastic ingestion on beta diversity in the gut. The expression of immune response genes was not altered in males. It did, however, cause some changes to alpha diversity metrics in both sexes and there were several differentially abundant taxa among treatments. These data suggest that microplastic ingestion has health effects, but these effects may be sex specific across certain species and they are likely not being solely driven by changes in gut microbial communities.

Microplastic Incubation:

We used 500-micron, blue, polyethylene beads (Cospheric, Inc., Santa Barbara, CA) for the incubation and ingestion study. First, a subsample of the new, clean microplastic beads were collected and stored at -80°C for 16S ribosomal RNA (rRNA) sequencing (‘pristine-MP’ sample). Then, the particles were secured in 16-micron nylon mesh bags attached to a PVC frame in Muskegon Lake, Michigan, USA at the Annis Water Resources Institute. The frame was secured at mid-depth in approximately 5 meters of water along the break wall from 14 June 2021 and incubated for a total of 56 days before final retrieval on 9 August 2021. Approximately every week during the incubation, the bags were inspected and gently rinsed and on 18 July 2021, after 34 days of incubation, we rotated the bags and examined them closely for wear, reinforcing seams as necessary. At final retrieval, the microplastics were removed gently from the incubation bags and combined in a composite container. We homogenized the weathered microplastics stored them in a 100 mL amber glass jar in Muskegon Lake surface water at 4°C for incorporation into the fish food for use in the ingestion study.

Addition to Food:

To prepare the fish food, individual 1.5 mm, slow-sinking, Finfish Starter pellets (Ziegler Bros Inc., Gardners, PA) were placed in wells of a standard 96-well PCR plate and moistened with laboratory grade deionized water. For treatment pellets, one microplastic bead was added to each well using fine-tipped forceps. For control pellets, no microplastic was added. We then pressed each softened pellet into the bottom of its well, embedding the microplastic particle into the treatment food pellets. As each plate was prepared, it was loosely covered by paper and placed in a dedicated refrigerator to dry overnight at 4°C. Pellets were divided into individual fish replicate meals in dedicated 5 mL tubes at the appropriate ratios for each treatment (outlined below). Pellets were prepared up to 4 days in advance and prepared meals stored at 4°C until use.

Ingestion Study Initiation:

The study was conducted in a flow-through system using Muskegon Lake surface water at the Mesocosm Lab of Annis Water Resources Institute in Muskegon, Michigan, USA. The study organisms were approximately five-month-old, adult fathead minnows hatched 16-20 February 2021 (Aquatic BioSystems, Inc.; Boulder, Colorado). The fish were received 21 July 2021 and acclimatized to the study system for 28 days before test initiation in 10-gallon glass aquaria in groups of approximately 24 fish with an equal distribution of males and females. Test chambers were made from bisected 10-gallon aquaria divided with plexiglass and opaque coverings between chambers. Each pair of chambers consisted of 1 male and 1 female of the same treatment. Chambers for each treatment group were randomly distributed throughout the test set up which consisted of three heavy-duty wire racks with two shelves each.

Fish were randomly distributed to individual test chambers four days prior to test initiation, on 7 August 2021 (day -4). Each treatment consisted of 14 individuals (7 males and 7 females). Initial measurements for standard length (SL), total length (TL), and mass were taken for each fish as they were distributed to their respective chamber. For the next three days leading up to the test initiation (days -3 to -1), pre-test fecal pellets were collected from the bottom of each chamber within 4 hours of being fed (‘before fecal’ samples). Each of these samples was gently centrifuged at 1000 rpm for approximately 30 seconds, the overlying water removed, and the fecal pellet frozen at -80°C for future DNA extraction for 16S rRNA sequencing of the before fecal microbial community. The exposure test began on 11 August 2021 (day 0). Each replicate was fed once daily at a rate of 3% by mass of the average body weight of all the test fish, which equaled an average of 16 pellets per day per fish. The test treatments were control, low, and high. Control replicates received 16 control pellets; low replicates received 4 microplastic-laden pellets and 12 control pellets; and high replicates received 16 microplastic-laden pellets each meal.

Chamber water level was maintained at approximately 10 liters with 2-4 daily water volume replacements maintained with adjustable flow valves on each chamber. Temperature was maintained at 23°C ± 2°C and dissolved oxygen (DO) was maintained at ≥ 6.0 mg/L. Temperature and DO were monitored in a subset of chambers daily while pH and nitrogen levels were monitored in the incoming water weekly. Feces, uneaten food, and other particulate was removed from each chamber by syphoning and disposed of daily before feeding and the entire system was scrubbed and flushed weekly to maintain water quality. No mortality was observed over the course of the study.

 Test Termination:

For the three days leading up to test termination (days 25-27), post-test fecal samples were collected from the bottom of each tank within four hours of the daily feedings (‘after fecal’ samples); samples were processed the same way as the pre-test fecal samples. The test was terminated on 8 September 2021 (day 28) with a random order for replicate processing. Each fish was measured for TL, SL, and mass in the same procedure as the pre-test measurements before being placed in a 500 mg/L solution of tricaine methanesulfonate (MS-222), adjusted to a pH between 7 and 8. After at least 5 minutes and complete anesthesia, each fish was taken from the solution, the caudle peduncle removed with a sterile scalpel and blood collected with a heparinized hematocrit tube. Each tube was plugged with clay and held vertically until centrifugation and measurement.

Each fish was then decapitated, and the liver was immediately removed, placed in  RNA-later (Invitrogen, Waltham, MA) and flash frozen in liquid nitrogen for future qPCR analysis of hepatic stress gene expression. The gastrointestinal tract from just posterior to the stomach to the anal pore was removed and the gut contents were expressed and placed into a 1.5 mL microcentrifuge tube (‘internal fecal’ samples) and the gut tissue was placed in another tube (‘gut tissue’ samples). Both the internal fecal and gut tissue samples were flash frozen with liquid nitrogen. All tools were flame sterilized and wiped down with RNA away (Thermo Scientific, Waltham, MA) between replicates. After each third of the test was taken down (14 replicates), all of samples were removed to the -80°C freezer and the hematocrit samples processed. Blood samples were centrifuged with a Micro-Hematocrit Centrifuge (Unico, Franksville, WI) at 11,500 rpm for 5 minutes. Readings were taken manually with calipers, measuring the length of the packed red blood cells and the whole blood.

 Sample Processing:

Gastrointestinal Microbial Community:

DNA was extracted from all 16S rRNA samples using PowerFecal Pro column extraction kits (Qiagen, Germantown, MA) according to the manufacturer protocol. A positive microbial community (ZymoBIOMICS Microbial Community Standards - cat. D6305, Zymo Research, Irvine, CA) and extraction blank were also extracted. Extracted samples were run on a 1% aragose gel for 60 minutes to verify genomic DNA was successfully extracted. Once all samples were extracted, the library was prepared according to the Illumina MiSeq 16S Metagenomic Sequencing Library Preparation protocol, with the amplification targeting the V4 region using the 515F and 806R primers (Caporaso et al., 2012). Samples were indexed using Nextera XT Index Primers (Illumina, Inc., San Diego, CA). After the final clean up step using the AMPure XP beads (Beckman Coulter, Brea, CA), each sample was quantified on a Qubit 3.0 Fluorometer using a dsDNA HS Assay kit (Invitrogen) to determine DNA concentration. Each sample was also run on a 2100 Bioanalyzer using a High Sensitivity DNA kit (Agilent, Santa Clara, CA) to determine target amplicon length and quality of the sample. Samples showing degradation were excluded from normalization and pooling. Acceptable samples were each normalized to 4 nM and pooled together. The sample library was denatured with freshly prepared 0.2N NaOH and diluted to 7 pM. PhiX was diluted to a concentration of 7 pM, denatured, and a 15% spike-in was included with the final sequencing library. Sequencing was performed on an Illumina MiSeq sequencer using a 2 x 250 bp format.

Gene Expression:

Liver RNA was extracted from male fathead minnows with a modified TRIzol protocol (Invitrogen). Each sample was thawed at room temperature until the liver tissue could be removed from the RNA later. The liver was then blotted dry and placed in a safe-lock microcentrifuge tube containing two 5 mm steel beads and placed in liquid nitrogen. The samples were then removed from the liquid nitrogen to allow the tissue to soften at room temperature for 3 minutes and homogenized in a TissueLyzer LT (Qiagen, Hilden, Germany) at 50 oscillations/s for 2 minutes. After phase separation with TRIzol/chloroform, and precipitation with isopropanol, the pellet was washed twice with 75% ethanol and resuspended in 100 μL of DEPC treated nuclease-free water.

The raw extracted RNA was quantified on a NanoDrop spectrophotometer (ThermoFisher, Waltham, MA) and the samples diluted to a target of 150 ng/μL at 100 μL, which was split into two 50 μL aliquots. One of these aliquots was cleaned using a Turbo DNA-free Kit (Invitrogen) and total RNA concentration, 260/280, and 260/230 ratios were assessed by NanoDrop after cleaning. Each sample was run on the 2100 Bioanalyzer using a Pico RNA kit (Agilent), to assess the quality of the extracted RNA, and quantified on a Qubit 3.0 Fluorimeter using an RNA HS/BR Assay Kit (Invitrogen) . Five samples from each treatment were selected for qPCR analysis based on quality as assessed through the bioanalyzer and these were diluted to 10 ng/μL.

Primers and probes were developed in house using the IDT PrimerQuest Tool (https://www.idtdna.com/pages/tools/primerquest) and sequences obtained from GenBank (Supplementary Table 1). Amplified products from each primer pair were sequenced on a 3500 Genetic Analyzer (Applied Biosystems, Waltham, MA) to verify correct target amplification before further validation. Validation included running a no reverse transcriptase control with each gene to verify the lack of DNA contamination remaining in the sample, testing the housekeeping gene across treatments to assess its constancy in expression, and calculating efficiency for each gene with a standard curve.

Hepatic gene expression was assessed with RT-qPCR on an Applied Biosystems OneStepPlus instrument. Gene targets were cyp1a, gst, and ncf2, which were normalized against elfa1 as the internal housekeeping gene. Relative gene expression was assessed using a ∆∆C_t protocol (Livak and Schmittgen, 2001; Schmittgen and Livak, 2008) on a OneStepPlus (Applied Biosystems) using Luna® Universal Probe One-Step RT-qPCR Kit (New England BioSystems). We used 2-μL of RNA template per well with each sample and control run in triplicate. Each reaction volume was 20 μL consisting of 10 μL of 2x Luna Universal Probe One-Step Reaction Mix, 1 μL of 20x Luna WarmStart RT Enzyme Mix, 0.4 uM final concentration of the  forward and reverse primers, 0.2 μM final concentration of the probe (FAM labeled with a Zen/Iowa BlackTM FQ double quencher), and the template RNA. The thermocycle for all validation and experimental plates for the gene expression analysis included an initial hold at 55°C for 10 min., followed by a denaturing step at 95°C for 1 min., then 40 cycles of 95°C for 10 sec., and 55°C for 1 min.

 Statistical Analysis:

Health Endpoints:

Change in mass was calculated as final mass (g) minus initial mass (g) for each replicate. Fulton’s Condition Factor (FCF) was calculated as the mass (g) divided by cubic length (cm). Change in FCF was calculated as final FCF minus initial FCF for each replicate. The hematocrit was calculated as the ratio of the red blood cells to the whole blood after centrifugation (Shiver & Grove, 2011). Exploratory data analysis showed interaction for all three of these endpoints between the sex and treatment factors and so growth, change in FCF, and hematocrit were analyzed for males and females separately. Shapiro-Wilk tests were used to test for normal distributions and Levene’s tests for equal variances. If all three treatments were normally distributed and had equal variances, they were compared by one-way ANOVA tests with TukeyHSD post-hoc. If any of the treatments were not normally distributed or had unequal variances, they were compared with Kruskal-Wallis tests with Wilcoxon post-hoc with Bonferroni adjustment. All analyses for health endpoints were conducted in R, version 4.2.0 (R Core Team, 2021).

Gastrointestinal Microbial Community

Sequence data were processed using the dada2 package, version 1.24.0, in R to generate amplicon sequence variants (ASVs) (Callahan et al., 2016). Briefly, sequences were filtered by quality scores, the sequencing error was ‘denoised’ by modeling error in the sequencing reads, chimeras were removed, and then the pairs merged. Sequences were aligned against the SILVA 16S reference database (Yilmaz et al., 2014). The resultant ASV and taxonomy tables were then analyzed in R using the phyloseq, version 1.40.0 (McMurdie and Holmes, 2013), vegan, version 2.6.2 (Oksanen et al., 2019), and DESeq2, version 1.36.0 (Love et al., 2014) packages. Phyloseq was used to compile taxonomy, ASV, and sample data into a single experiment-level object. Exploratory data analysis on all samples included assessing family level composition with stacked bar charts and unconstrained ordinations with principal coordinate analysis (PCoA).

Samples were split into internal fecal, gut tissue, before fecal, and after fecal sample types and then each of those was subset into males and females for independent analysis. To reiterate, ‘internal fecal’ samples consisted of the fecal material removed from the gut at test termination, ‘gut tissue’ samples consisted of the remaining gastrointestinal tract after feces was removed at test termination, ‘before fecal’ samples consisted of fecal material sampled from the bottom of the tank after initial acclimation but prior to the start of the ingestion study, and ‘after fecal’ samples consisted of fecal material collected from the bottom of the tank at the end of the ingestion study. We used the control group of the internal fecal samples for males and females as a representative subset of samples to assess sex differences.

We assessed the community composition at the beginning of the experiment using the before fecal samples in each sex and expected no differences between treatments since all individuals were on the same standardized diet prior to the beginning of the study.  To investigate whether the gut microbial community changed over the course of the ingestion period, we compared the before fecal and after fecal samples within each treatment. We used the internal fecal, gut tissue, and after fecal samples to assess the differences between treatments after the microplastic exposure period.

For diversity metrics, each group was rarefied to the lowest sample read depth in each group (Supplemental Table 2). For Chao1 richness estimates and Shannon diversity estimates, Shapiro-Wilk tests were used to test for normal distributions and Levene’s tests for equal variances. If all levels were normally distributed and had equal variances, they were compared by one-way ANOVA tests with TukeyHSD post-hoc. If any of the treatments were not normally distributed or had unequal variances, they were compared with Kruskal-Wallis tests with Wilcoxon post-hoc with Bonferroni adjustment or, in the case of the fecal control group, the Wilcoxon rank sum test. Chao1 and Shannon estimates were also calculated and compared for the before and after samples together to assess changes over time in each treatment. As these are paired data, they were tested with either paired T-tests or Wilcoxon signed-rank tests. Beta diversity was explored by generating PCoA ordinations based on a Bray-Curtis distance matrices within each sample type. Permutational multivariate analysis of variance (PERMANOVA) tests were performed with the adonis function in the vegan package for each sex in each sample type. As the before and after fecal samples assessed together violated the assumption of independence, these were not tested by PERMANOVA. Differential abundance was assessed on non-rarefied count data using DESeq2 which assumes a negative binomial distribution of the abundance count data. Data was independently filtered to include only those taxa present at least three times and in at least 10% of the samples. The DESeq2 results were then subject to log fold change shrinking with the lfcShrink function to account for effect sizes. Taxa were considered biologically significant in their differential abundances if they had a log2 fold change value greater than 0.25 and a p-value less than 0.05. This threshold translates to approximately +/- 20% in abundance and allowed us to focus on those taxa with the greatest differences among treatment types; those which may have had impacts on the microbial community function. Raw sequencing data has been deposited on NCBI’s Sequence Read Archive (SRA) with accession number SUB14107689.

Gene Expression

The coefficient of variation (CV) was calculated for the housekeeping gene, ef1a, using transformed C_t values () for all study replicates. For each target gene (cyp1a, gsta, and ncf2) the coefficient of variation was also calculated for the controls using the transformed ∆C_t values (), as calculated in Equation 1 in manuscript.

As there was high variation in the controls for gsta in particular, we used the geometric mean to calculate the ∆∆C_t values. Gene expression levels for each of the three target genes were normalized against the housekeeping gene efla1 using the ∆∆C_t method, Equation 2 (in manuscript)

These values were then transformed into gene fold, Equation 3 (in manuscript)

These gene fold values were log-transformed then compared across treatment groups with one-way ANOVA for each gene.