diff --git a/DESCRIPTION b/DESCRIPTION
index 2b9e5e3cd..f4280d4ef 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,5 +1,5 @@
 Package: AMR
-Version: 3.1.0.9018
+Version: 3.1.0.9019
 Date: 2026-11-09
 Title: Antimicrobial Resistance Data Analysis
 Description: Functions to simplify and standardise antimicrobial resistance (AMR)
diff --git a/NEWS.md b/NEWS.md
index d47323807..384d878e6 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,4 @@
-# AMR 3.1.0.9018
+# AMR 3.1.0.9019
 
 This version fixes a critical bug in WISCA that biased all coverage estimates of `wisca()` and `antibiogram(..., wisca = TRUE)` downwards since AMR 3.0.0. Users are strongly advised to re-run their WISCA analyses.
 
@@ -29,6 +29,7 @@ It also adds a `site` argument to `as.sir()`, stores EUCAST MIC breakpoints exac
 * **Taxonomic names in plots**
   * New function `labels_taxonomy()` to print taxonomic names in axis and legend labels in italics, which can directly be used as labeller in `ggplot2` (e.g., `scale_x_discrete(labels = labels_taxonomy)`) and in base R plots
   * `italicise_taxonomy()` gained `type = "plotmath"`, which returns an expression vector for titles, axes and legends in plots
+  * All plotting functions now print taxonomic names in italics using `labels_taxonomy()`: the facets and titles of `autoplot()` and `plot()` for antibiograms, the legend of `wisca_plot()` (which previously italicised all its entries, including e.g. "CoNS"), the axes, legends and facets of `ggplot_sir()` and `facet_sir()`, the labels and legend of `ggplot_pca()`, and the subtitles of `plot()` and `autoplot()` for MIC and disk diffusion values (which previously italicised the whole name, e.g. "Coagulase-negative Staphylococcus (CoNS)")
 
 ### Changed
 * SIR interpretation (`as.sir()`)
@@ -43,6 +44,8 @@ It also adds a `site` argument to `as.sir()`, stores EUCAST MIC breakpoints exac
   * Fixed punctuation being italicised along with taxonomic names, e.g., `"(S. aureus)"` now gives `"(*S. aureus*)"` instead of `"(S. *aureus)*"`
   * `NA` values now return `NA` instead of `"NA"`
   * Taxonomic names are now looked up once per session, making the function considerably faster for long vectors
+* `ggplot_pca()`
+  * Fixed an error when setting `groups` together with `labels`, and setting `groups` without `labels` now colours the points by these groups instead of using them as labels, as documented
 * *Morganella* is interpreted as reported, with the breakpoints of Enterobacterales or Morganellaceae and not of *M. morganii*, which is now tested (#330)
 * WISCA (`wisca()` and `antibiogram(..., wisca = TRUE)`)
   * Fixed a critical bug that systematically underestimated the coverage of every WISCA since AMR 3.0.0. Pathogens absent from a syndromic group, and pathogens without any test result for a regimen, were still included in the model with only prior information: each added a prior weight of 1 to the pathogen incidence, and untested pathogens kept their full incidence weight while counting as about 50% susceptible. The bias was largest for stratified analyses (e.g., up to 14 percentage points for meropenem in random strata of `example_isolates`). Now, only pathogens that occur in a group and have at least one test result for the regimen are included. Please re-run all WISCA analyses (#337)
diff --git a/R/antibiogram.R b/R/antibiogram.R
index bc1aac90d..188ca4609 100755
--- a/R/antibiogram.R
+++ b/R/antibiogram.R
@@ -1404,7 +1404,7 @@ plot.antibiogram <- function(x, ...) {
       names.arg = df_sub$ab,
       col = "#aaaaaa",
       beside = TRUE,
-      main = mo,
+      main = labels_taxonomy(mo),
       legend = NULL
     )
 
@@ -1568,7 +1568,7 @@ autoplot.antibiogram <- function(object,
   }
 
   if (has_facet) {
-    out <- out + ggplot2::facet_wrap("mo")
+    out <- out + ggplot2::facet_wrap("mo", labeller = labeller_taxonomy)
   }
 
   out
@@ -1733,6 +1733,7 @@ plot_wisca_susceptibility_incidence <- function(wisca_model, sep) {
       colour = translate_AMR("Pathogen", language = get_AMR_locale()),
       caption = paste(attributes(wisca_model)$simulations, "Monte Carlo simulations")
     ) +
+    ggplot2::scale_colour_discrete(labels = labels_taxonomy) +
     ggplot2::guides(
       colour = ggplot2::guide_legend(
         override.aes = list(alpha = 1, size = 3)
@@ -1741,7 +1742,6 @@ plot_wisca_susceptibility_incidence <- function(wisca_model, sep) {
     ggplot2::theme(
       legend.position = "right",
       legend.key.spacing.y = ggplot2::unit(0.25, "lines"),
-      legend.text = ggplot2::element_text(face = "italic"),
       plot.title = ggplot2::element_text(size = 12, face = "bold")
     )
 }
diff --git a/R/ggplot_pca.R b/R/ggplot_pca.R
index 89225a6f7..2622a79db 100755
--- a/R/ggplot_pca.R
+++ b/R/ggplot_pca.R
@@ -151,6 +151,7 @@ ggplot_pca <- function(x,
     pca_model = x,
     groups = groups,
     groups_missing = missing(groups),
+    group_name = if (missing(groups)) NULL else paste(deparse(substitute(groups)), collapse = ""),
     labels = labels,
     labels_missing = missing(labels),
     choices = choices,
@@ -183,8 +184,19 @@ ggplot_pca <- function(x,
   )
 
   # Score Labels
+  labels_parse <- FALSE
   if (!is.null(labels)) {
     df.u$labels <- labels
+    if (contains_taxonomy(labels)) {
+      # italicise taxonomic names, as plotmath text to be parsed by geom_text()
+      labels_parse <- TRUE
+      df.u$labels <- vapply(
+        FUN.VALUE = character(1),
+        as.list(labels_taxonomy(labels)),
+        function(e) if (is.character(e) && is.na(e)) NA_character_ else paste(deparse(e, width.cutoff = 500L), collapse = ""),
+        USE.NAMES = FALSE
+      )
+    }
   }
 
   # Grouping variable
@@ -213,7 +225,8 @@ ggplot_pca <- function(x,
       ) +
         ggplot2::geom_text(ggplot2::aes(label = labels, colour = groups),
           nudge_y = -0.05,
-          size = labels_textsize
+          size = labels_textsize,
+          parse = labels_parse
         ) +
         ggplot2::labs(colour = group_name)
     } else {
@@ -223,7 +236,8 @@ ggplot_pca <- function(x,
       ) +
         ggplot2::geom_text(ggplot2::aes(label = labels),
           nudge_y = -0.05,
-          size = labels_textsize
+          size = labels_textsize,
+          parse = labels_parse
         )
     }
   } else {
@@ -241,6 +255,11 @@ ggplot_pca <- function(x,
     }
   }
 
+  # italicise taxonomic names in the legend
+  if (!is.null(df.u$groups) && contains_taxonomy(df.u$groups)) {
+    g <- g + ggplot2::scale_colour_discrete(labels = labels_taxonomy)
+  }
+
   # Overlay a concentration ellipse if there are groups
   if (!is.null(df.u$groups) && !is.null(ell) && isTRUE(ellipse)) {
     g <- g + ggplot2::geom_path(
@@ -308,6 +327,7 @@ ggplot_pca <- function(x,
 pca_calculations <- function(pca_model,
                              groups = NULL,
                              groups_missing = TRUE,
+                             group_name = NULL,
                              labels = NULL,
                              labels_missing = TRUE,
                              choices = 1:2,
@@ -329,8 +349,8 @@ pca_calculations <- function(pca_model,
       error = function(e) NULL
     )
   }
-  if (!is.null(groups) && is.null(labels)) {
-    # turn them around
+  if (groups_missing && !is.null(groups) && is.null(labels)) {
+    # only one non-numeric column: use it for labels instead of groups
     labels <- groups
     groups <- NULL
     group_name <- NULL
diff --git a/R/ggplot_sir.R b/R/ggplot_sir.R
index 3bb5fc1f7..153af875f 100755
--- a/R/ggplot_sir.R
+++ b/R/ggplot_sir.R
@@ -270,6 +270,14 @@ ggplot_sir <- function(data,
     )
   }
 
+  # italicise taxonomic names on the x axis and in the legend, e.g. with x = "mo" or fill = "mo"
+  if (x %in% colnames(data) && contains_taxonomy(data[[x]])) {
+    p <- p + ggplot2::scale_x_discrete(labels = labels_taxonomy)
+  }
+  if (fill != "interpretation" && fill %in% colnames(data) && contains_taxonomy(data[[fill]])) {
+    p <- p + ggplot2::scale_fill_discrete(labels = labels_taxonomy)
+  }
+
   if (!is.null(facet)) {
     # pass the value, as facet_sir() takes an unquoted name as column name
     p <- p + do.call(facet_sir, list(facet = facet, nrow = nrow))
diff --git a/R/italicise_taxonomy.R b/R/italicise_taxonomy.R
index 2abca5693..d5ea9baf8 100755
--- a/R/italicise_taxonomy.R
+++ b/R/italicise_taxonomy.R
@@ -107,6 +107,12 @@ italicise_taxonomy <- function(string, type = c("markdown", "ansi", "html", "plo
   )
 }
 
+contains_taxonomy <- function(x) {
+  # whether any value holds a taxonomic name that would be italicised
+  x <- unique(as.character(x))
+  any(vapply(taxonomy_pieces(x[!is.na(x)]), function(p) any(p$italic), logical(1)))
+}
+
 get_taxonomy_names <- function() {
   # cached for the session, reset by reset_mo_cache() when the lookup table changes
   if (is.null(AMR_env$taxonomy_names)) {
diff --git a/R/plotting.R b/R/plotting.R
index a8e4bb8ce..a8fec12bf 100755
--- a/R/plotting.R
+++ b/R/plotting.R
@@ -80,6 +80,8 @@
 #'
 #' The function [labels_taxonomy()] prints taxonomic names in axis and legend labels in italics, as is the convention in microbiology, e.g., *E. coli* and *Klebsiella pneumoniae*. It can directly be used as labeller in `ggplot2`, e.g., `scale_x_discrete(labels = labels_taxonomy)`, and in base \R, e.g., `axis(1, at = 1:3, labels = labels_taxonomy(names))`. Other words, such as "Other" or "CoNS", stay upright. It uses [italicise_taxonomy()] with `type = "plotmath"`, which also works for titles. It returns an [expression] vector of the same length as `x`; other arguments in `...` are ignored.
 #'
+#' The plotting functions of this package apply this automatically wherever they show taxonomic names, e.g., in the facets and titles of [antibiogram()] plots, in the legend of [wisca_plot()], on the axes, legends and facets of [ggplot_sir()] and [facet_sir()], in [ggplot_pca()], and in the subtitles of MIC and disk diffusion plots.
+#'
 #' The interpretation of "I" will be named "Susceptible, incr. exp." (susceptible, increased exposure) for all EUCAST guidelines since 2019, and will be named "Intermediate" in all other cases.
 #'
 #' For interpreting MIC values as well as disk diffusion diameters, the default guideline is `r AMR::clinical_breakpoints$guideline[1]`, unless the package option [`AMR_guideline`][AMR-options] is set. See [as.sir()] for more information.
@@ -1288,7 +1290,7 @@ plot_colours_subtitle_guideline <- function(x, mo, ab, guideline, colours_SIR, f
     sir <- as.character(sir)
     cols <- unname(plot_sir_colour_values(colours_SIR)[sir])
     cols[is.na(cols)] <- "#BEBEBE"
-    sub <- bquote(.(abname) ~ "-" ~ italic(.(moname)) ~ .(guideline_txt))
+    sub <- bquote(.(abname) ~ "-" ~ .(italicise_taxonomy(moname, type = "plotmath")[[1]]) ~ .(guideline_txt))
   } else {
     sir <- rep(NA_character_, length(x))
     cols <- "#BEBEBE"
@@ -1319,7 +1321,7 @@ facet_sir <- function(facet = c("interpretation", "antibiotic"), nrow = NULL) {
     facet <- "antibiotic"
   }
 
-  ggplot2::facet_wrap(facets = facet, scales = "free_x", nrow = nrow)
+  ggplot2::facet_wrap(facets = facet, scales = "free_x", nrow = nrow, labeller = labeller_taxonomy)
 }
 
 #' @rdname plot
@@ -1437,6 +1439,17 @@ labels_taxonomy <- function(x, ...) {
   italicise_taxonomy(as.character(x), type = "plotmath")
 }
 
+# ggplot2 facet labeller that italicises taxonomic names in the strips
+labeller_taxonomy <- structure(
+  function(labels) {
+    lapply(labels, function(values) {
+      expr <- labels_taxonomy(values)
+      lapply(seq_along(expr), function(i) expr[i])
+    })
+  },
+  class = c("function", "labeller")
+)
+
 #' @rdname ggplot_sir
 #' @export
 labels_sir_count <- function(position = NULL,
diff --git a/man/plot.Rd b/man/plot.Rd
index bf6067202..aff08c050 100644
--- a/man/plot.Rd
+++ b/man/plot.Rd
@@ -238,6 +238,8 @@ This package contains more functions that extend the \code{ggplot2} package, to
 
 The function \code{\link[=labels_taxonomy]{labels_taxonomy()}} prints taxonomic names in axis and legend labels in italics, as is the convention in microbiology, e.g., \emph{E. coli} and \emph{Klebsiella pneumoniae}. It can directly be used as labeller in \code{ggplot2}, e.g., \code{scale_x_discrete(labels = labels_taxonomy)}, and in base \R, e.g., \code{axis(1, at = 1:3, labels = labels_taxonomy(names))}. Other words, such as "Other" or "CoNS", stay upright. It uses \code{\link[=italicise_taxonomy]{italicise_taxonomy()}} with \code{type = "plotmath"}, which also works for titles. It returns an \link{expression} vector of the same length as \code{x}; other arguments in \code{...} are ignored.
 
+The plotting functions of this package apply this automatically wherever they show taxonomic names, e.g., in the facets and titles of \code{\link[=antibiogram]{antibiogram()}} plots, in the legend of \code{\link[=wisca_plot]{wisca_plot()}}, on the axes, legends and facets of \code{\link[=ggplot_sir]{ggplot_sir()}} and \code{\link[=facet_sir]{facet_sir()}}, in \code{\link[=ggplot_pca]{ggplot_pca()}}, and in the subtitles of MIC and disk diffusion plots.
+
 The interpretation of "I" will be named "Susceptible, incr. exp." (susceptible, increased exposure) for all EUCAST guidelines since 2019, and will be named "Intermediate" in all other cases.
 
 For interpreting MIC values as well as disk diffusion diameters, the default guideline is EUCAST 2026, unless the package option \code{\link[=AMR-options]{AMR_guideline}} is set. See \code{\link[=as.sir]{as.sir()}} for more information.
diff --git a/tests/testthat/test-italicise_taxonomy.R b/tests/testthat/test-italicise_taxonomy.R
index 8c44e89b2..e4169dba5 100644
--- a/tests/testthat/test-italicise_taxonomy.R
+++ b/tests/testthat/test-italicise_taxonomy.R
@@ -73,4 +73,12 @@ test_that("test-italicise_taxonomy.R", {
     expression(italic("E. coli"), "Other")
   )
   expect_identical(labels_taxonomy(character(0)), expression())
+
+  # used internally by the plotting functions
+  expect_true(AMR:::contains_taxonomy(c("E. coli", "Other", NA)))
+  expect_false(AMR:::contains_taxonomy(c("Other", "CoNS", NA)))
+  expect_identical(
+    AMR:::labeller_taxonomy(list(mo = c("E. coli", "CoNS"))),
+    list(mo = list(expression(italic("E. coli")), expression("CoNS")))
+  )
 })
diff --git a/tests/testthat/test-pca.R b/tests/testthat/test-pca.R
index 0ca055de2..4be80f0f9 100644
--- a/tests/testthat/test-pca.R
+++ b/tests/testthat/test-pca.R
@@ -60,6 +60,13 @@ test_that("test-pca.R", {
   if (AMR:::pkg_is_available("ggplot2")) {
     ggplot_pca(pca_model, ellipse = TRUE)
     ggplot_pca(pca_model, arrows_textangled = FALSE)
+    # explicit groups, with and without labels
+    my_groups <- c("Gram-positive", "Gram-negative", "Gram-negative")
+    p <- ggplot_pca(pca_model, groups = my_groups)
+    expect_identical(p$data$groups, my_groups)
+    expect_identical(p$labels$colour, "my_groups")
+    p <- ggplot_pca(pca_model, labels = resistance_data$genus, groups = my_groups)
+    expect_identical(p$data$groups, my_groups)
   }
 
   if (AMR:::pkg_is_available("dplyr", min_version = "1.0.0", also_load = TRUE)) {