Yes, you can do such a thing in JuMP.
If you use the so-called direct mode, then incremental changes to the model are passed directly to the solver.
If your solver supports incremental changes and re-solve, then it will do so.
For instance, if you are solving an LP, adding a constraint, and re-solving, then the second solve will (likely) be a lot faster than the first.
All mainstream solvers I know of, e.g., CPLEX, Gurobi, GLPK, HiGHS, etc., support this.
BTW: nowadays, I would recommend using HiGHS instead of GLPK. It's much much faster and more numerically robust.
Here is a small example based on a knapsack problem
# Generate instance data
n = 8 # number of items
w = rand(n) # item weights
c = rand(n) # item values
W = sum(w) / 2 # total knapsack capacity
# Build the optimization model
# (we only give the LP relaxation to GLPK)
model = JuMP.Model(GLPK.Optimizer)
# Create x variables in [0, 1]
# they are not marked as binary to do branching manually later
@constraint(model, bounds, 0 .<= x .<= 1)
@constraint(model, capacity, w'x <= W) # knapsack capacity
@objective(model, Max, c'x) # maximize value of selected items
# Main loop
z_root = objective_value(model)
x_ = value.(x)
# branch on first fractional variable
# (here we assume that `x1` is fractional)
# Left (x1 = 0) node
z0 = objective_value(model)
# Right (x1 = 1) node
z1 = objective_value(model)
If you have more technical questions, I encourage you to post them on the Julia discourse forum on mathematical optimization. You'll get more traction there.