It seems that main issue was the initial guess.
The algorithm was not able to reach an equilibrium point from the provided initial guess.

I'm no expert, but with minor tweaks I managed to get the solutions you mentioned.

The code:

using Revise, Parameters, Plots
using BifurcationKit
const BK = BifurcationKit

# vector field of the problem
function COm(u, p)
    @unpack r,K,a,h,eps,mu = p
    x, y = u
    out = similar(u)
    out[1] = r*x*(1.0 - x/K) - a*x*y/(1 + a*h*x)
    out[2] = eps*a*x*y/(1 + a*h*x) - mu*y
    out
end

####
# integrate ODE for sanity check
# a=0.1 stable
# a=0.3 oscillatory
using DifferentialEquations
z0 = [0.05, 0.1]
alg_ode = Rodas5()
# stable
@reset par_com.a = 0.1
prob_de = ODEProblem(COm, z0, (0.0, 300.0), par_com)
sol_ode = solve(prob_de, alg_ode);
plot(sol_ode)
# keep equilibrium point
zEq = sol_ode.u[end]

# oscillatory
@reset par_com.a = 0.3
prob_de = ODEProblem(COm, z0, (0.0, 300.0), par_com)
sol_ode = solve(prob_de, alg_ode);
plot!(sol_ode)


# parameters used in the model
par_com = (r = 1.0, K = 10.0, a = 0.1, h = 0.5, eps = 0.5, mu = 0.2)
# record variables for plotting
recordCO(x, p) = (x = x[1], y = x[2])

# initial condition
# z0 = [0.05, 0.1] ## fail to find equilibrium
# z0 = [5.0, 6.3] ## guess, from numerical data
z0 = zEq # initiate from numerical data

# Bifurcation Problem
prob = BifurcationProblem(COm, z0, par_com, (@optic _.a); record_from_solution = recordCO)

# continuation parameters
opts_br = ContinuationPar(p_min = 0.01, p_max = 0.5, dsmin=0.01, ds = 0.1, dsmax = 1.0)

# compute the branch of solutions
br = continuation(prob, PALC(), opts_br; plot = true, verbosity = 0, normC = norminf,bothside = true)
# plot the branch
scene = plot(br)

## inspect transient before and after the special points
pt = 3 # point of interest
prob_de = ODEProblem(COm, br.specialpoint[pt].x, (0.0, 1000.0), 
    @set par_com.a = br.specialpoint[pt].param - 0.01)
sol_ode = solve(prob_de, alg_ode);
plot(sol_ode)
prob_de = ODEProblem(COm, br.specialpoint[pt].x, (0.0, 1000.0), 
    @set par_com.a = br.specialpoint[pt].param + 0.01)
sol_ode = solve(prob_de, alg_ode);
plot!(sol_ode)

The continuation result:

julia> br
 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 272
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter a starts at 0.01, ends at 0.5
 ├─ Algo: PALC
 └─ Special points:

- #  1, endpoint at a ≈ +0.01000000,                                                                     step =   0
- #  2,       bp at a ≈ +0.04983159 ∈ (+0.04983159, +0.05006280), |δp|=2e-04, [converged], δ = (-1,  0), step = 252
- #  3,     hopf at a ≈ +0.30509470 ∈ (+0.29935407, +0.30509470), |δp|=6e-03, [converged], δ = ( 2,  2), step = 269
- #  4, endpoint at a ≈ +0.50000000,  

I was going to add the images, but 10 reputation points are needed. It is my first post in here...