You can easily do it by modifying the plotPCA() function. My plotPCA is a modified one, which plots the sampleName along with the bubbles.

At Line 24 to 29 >>

  d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], group = group, 
                  intgroup.df, name = colData(rld)[,1])
  if (returnData) {
    attr(d, "percentVar") <- percentVar[1:2]
    return(d)
  }

Change PC1 and PC2 to pca$x[, 2] and pca$x[, 3] respectively. Accordingly, the percentVar[1:2] to percentVar[2:3]

From what I can tell there is no way to return PC2/PC3 from the plotPCA function. Instead what you could do is just use the code from within the plotPCA function to return the full results of the principal component analysis:

  # calculate the variance for each gene
  rv <- rowVars(assay(object))

  # select the ntop genes by variance
  select <- order(rv, decreasing=TRUE)[seq_len(min(ntop, length(rv)))]

  # perform a PCA on the data in assay(x) for the selected genes
  pca <- prcomp(t(assay(object)[select,]))

  # the contribution to the total variance for each component
  percentVar <- pca$sdev^2 / sum( pca$sdev^2 )