Change legend and shape in ggbiplot pca-CodePudding

can you please help me with my pca? I would like to change the shapes in that way that each species has a different color and all 2-3 organisms for each species have 2-3 different symbols. It should look like this:
DATA
results <- structure(list(organism = c(
  "Mt1", "Mt1", "Mt1", "Mt2", "Mt2",
  "Mt2", "Al1", "Al1", "Al1", "Al2", "Al2", "Al2", "Bg1", "Bg1",
  "Bg1", "Bg2", "Bg2", "Bg2", "Bg3", "Bg3", "Bg3", "Cs1", "Cs1",
  "Cs1", "Cs2", "Cs2", "Cs2", "Cs3", "Cs3", "Cs3", "Df1", "Df1",
  "Df1", "Df2", "Df2", "Df2", "Df3", "Df3", "Df3"
), species = c(
  "Mt",
  "Mt", "Mt", "Mt", "Mt", "Mt", "Al", "Al", "Al", "Al", "Al", "Al",
  "Bg", "Bg", "Bg", "Bg", "Bg", "Bg", "Bg", "Bg", "Bg", "Cs", "Cs",
  "Cs", "Cs", "Cs", "Cs", "Cs", "Cs", "Cs", "Df", "Df", "Df", "Df",
  "Df", "Df", "Df", "Df", "Df"
), lactones = c(
  23435.167, 23303.111,
  22340.556, 51519.222, 48824.5, 50427.278, 25260.222, 23177.556,
  27903, 5993.722, 7307.389, 9419.167, 58216.944, 69860.222, 72809.333,
  51584.611, 57738.056, 53356.056, 41983.389, 46930.722, 49487.611,
  7155.056, 8153.611, 7445.722, 10771.556, 12239.556, 13788.611,
  12082.111, 14331.556, 14123.056, 26906.778, 20689.111, 22872.722,
  18838.222, 18218.667, 18545.389, 19924.111, 22504.056, 16907.278
), cyanides = c(
  166.4, 168.9, 176.6, 175.9, 166.5, 165.4, 162,
  161.5, 156.2, 180.4, 169.7, 147.5, 132.4, 147.4, 145.8, 142.9,
  141.2, 142.7, 130.8, 148.3, 139.4, 161.6, 151, 168.6, 126.4,
  142.6, 136.1, 152, 143.3, 158.2, 156.2, 163.5, 154.1, 159.2,
  155.8, 156.2, 156.4, 158, 151.5
), alcohols = c(
  137653.9, 153511,
  150719.9, 173401.1, 171614.4, 168865.1, 211737.4, 199886.7, 240088.4,
  289334.9, 275631.1, 277924.5, 92275.3, 105654.9, 111731.3, 105548.2,
  117728.9, 110260.2, 103799.4, 113944.6, 121976.5, 221538.8, 179823,
  176978, 144314.1, 142913.7, 131506.7, 171730, 141748.7, 150303,
  310203.9, 214322.9, 197572.9, 125167.6, 127077.2, 154400.4, 199472.6,
  196343, 185052.9
), ethers = c(
  4040.1, 4529.1, 3255.8, 3890.1,
  3694.1, 3693.2, 9563.9, 10403.3, 11897.1, 6673.3, 8333.4, 9622.2,
  4099.5, 4080.6, 4014.6, 6450.1, 6332.4, 6506.2, 4781.8, 5151.6,
  5318.3, 8356, 7961.2, 8196.5, 8634.6, 9471.9, 9390.4, 6259.6,
  7532.8, 7755.8, 5505.5, 5315.6, 4627.7, 12372.9, 11182, 10543.6,
  4452.3, 3994.1, 4084.6
), acids = c(
  1131.52, 1148.52, 1200.88,
  1196.12, 1132.2, 1124.72, 1101.6, 1098.2, 1062.16, 1226.72, 1153.96,
  1003, 900.32, 1002.32, 991.44, 971.72, 960.16, 970.36, 889.44,
  1008.44, 947.92, 1098.88, 1026.8, 1146.48, 859.52, 969.68, 925.48,
  1033.6, 974.44, 1075.76, 1062.16, 1111.8, 1047.88, 1082.56, 1059.44,
  1062.16, 1063.52, 1074.4, 1030.2
)), class = "data.frame", row.names = c(
  "1",
  "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",
  "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24",
  "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35",
  "36", "37", "38", "39"
))